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AN ITERATIVE TRANSFORMATION PROCEDURE FOR NUMERICAL SOLUTION OF FLUTTER 

AND SIMILAR CHARACTERISTIC -VALUE PROBLEMS 1 

By Myeon L. Gossahd 


SUMMARY 

An iterative transformation procedure suggested by H. 
Wielandt for numerical solution of fluff, er and similar charac- 
teristic-value problems is presented. Application of this pro- 
cedure to ordinary natural-vibration problems and to flutter 
problems is shown by numerical examples. Comparisons of 
computed results with experimental values and with results 
obtained by other methods of analysis are made. 

INTRODUCTION 

Existing methods of flutter analysis include the 
representative-section method, generalized-coordinate meth- 
ods, matrix methods, and operational methods. The present 
report presents an iteration procedure for analysis of 
flutter and similar characteristic-value problems. 

For ordinary natural-vibration problems, iterative pro- 
cedures of the Stodola type (references 1 and 2) are suitable 
for finding the fundamental and higher-order natural modes 
and frequencies. The higher-order solutions are obtained 
through use of the orthogonality relations that exist among 
the natural modes. 

Orthogonality relations analogous to those that exist in 
ordinary vibration problems can be found in the flutter 
problem only by introduction of the so-called “adjoint” 
problem. (This additional step is unnecessary in ordinary 
vibration problems by virtue of the fact that they are 
“self-adjoint.”) Wielandt has suggested an iterative trans- 
formation procedure (reference 3) which is well-suited to 
the flutter problem in that it avoids the need of orthogo- 
nality relations and hence does not require consideration of 
the adjoint problem. The iterative transformation pro- 
cedure can also be applied to ordinary natural-vibration 
problems with less labor than is generally required in the 
extended Stodola procedure. 

Because both the original and translated works of 
Wielandt are difficult to follow, an explanation of the idea 
of the iterative transformation procedure is given in the 
present report and the application of the procedure to 
ordinary natural-vibration problems and to flexure-torsion 
flutter problems is shown in numerical examples. Com- 
parisons of computed results with experimental values and 
with results obtained by other methods of analysis are 
also made. 


SYMBOLS 

El flexural stiffness 

GJ torsional stiffness 

x span wise coordinate with origin at root of wing 

y complex representation of amplitudes and phases 

of translation of elastic axis in harmonic motion 
< p complex representation of amplitudes and phases 

of rotation about elastic axis in harmonic motion 
w coupled mode (?/,<£) 

P v complex coefficients of y which, when multiplied 

by y, give complex representation of amplitudes 
and phases of aerodynamic and inertia forces 
associated with translational component of 
harmonic motion 

P$ complex coefficients of 4> which, when multiplied 

by <£, give complex representation of amplitudes 
and phases of aerodynamic and inertia forces 
associated with rotational component of har- 
monic motion 

Q v complex coefficients of y which, when multiplied 

by y, give complex representation of amplitudes 
and phases of aerodynamic and inertia torques 
associated with translational component of har- 
monic motion 

complex coefficients of <j> which, when multiplied 
by <p, give complex representation of amplitudes 
and phases of aerodynamic and inertia torques 
associated with rotational component of har- 
monic motion 

g v ,g,j, structural-damping coefficients associated with 

flexure and torsion, respectively (see appendix B) 
g a coefficient of artificial damping (may be either 

positive or negative) 
k reduced frequency (bu/v) 

co frequency of harmonic motion 

C characteristic value 

b length of semichord of wing 

L length of cantilever wing from root to tip 

H mass ratio (y/irph 2 ) 

v velocity of air relative to wing 

7 distributed mass of wing per unit length of span 

p mass density of air 


(rP) 


1 Supersedes NACA TN 2346, “An Iterative Transformation Procedure for Numerical Solution of Flutter and Similar Characteristic-Value Problems” by Myron L. Gossard, 1951. 
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a distance between midchord axis and elastic axis 

in terms of local semichord, positive when 
elastic axis is behind midchord axis 
u distance between elastic axis and gravity axis of dis- 

tributed mass of wing in terms of local semichord, 
positive when gravity axis is behind elastic axis 
r radius of gyration of distributed mass of wing 

about elastic axis in terms of local semichord 
F,G transcendental functions of k (see reference 4) 
t time 

F mn eigenvalue factor ^^—1^ 

R ratio of complex constants 

X length; in numerical solutions, distance between 

specific adjacent stations of wing 
p applied force 

g applied torque 

V shear 

M bending moment 


a curvature 

/3 slope of elastic axis 

T twisting moment 

6 angle of twist 

Subscripts : 

1,2,3, . . . 
a2,<z3,a4, . . . 
b 

A, BA • ■ • 

R 
I 
0 

bl,ba2,ba3, . . . 


true modes or eigenvalues 
transformed modes 
intermediate derived mode 
stations 
real 

imaginary 
reference value 
sweeping functions 


Superscripts : 

(1),(2),(3), . . . cycles of iteration 

A bar over a symbol indicates a concentrated quantity 
instead of a distributed quantity. 

A prime is used to denote division by w 2 . 


ITERATIVE TRANSFORMATION METHOD OF SOLUTION 

GENERAL FEATURES OF METHOD 


The principle of the iterative transformation procedure is 
similar in form to that of the standard iteration procedure 
for solving characteristic-value problems. Both procedures 
require the determination of the solutions in the order of the 
magnitudes of the eigenvalues, beginning with the funda- 
mental. Both procedures require assumptions of modes, 
integrations which generally must be done numerically, and 
sweeping operations for higher-order-mode determinations. 
The distinguishing features of the iterative transformation 
procedure occur in the determination of solutions higher than 
the fundamental and are as follows: (1) The immediate aim 
is to determine not the true wth mode, as in the standard 
iteration procedure, but a particular linear combination 
composed of all modes from the fundamental to the ■rath. 
This linear combination is referred to as the transformed nth 


mode. The transformed nth mode can be made to have 
nodal (zero) points at specified stations of the wing; such a 
feature is highly desirable in numerical work. (2) The 
sweeping operations, which consist of subtractions of lower- 
order-mode shapes from the function obtained by integrating 
the assumed mode, do not employ the orthogonality relations 
as in the standard iteration method but make use of forcing 
functions that, in numerical work, greatly simplify the 
sweeping operations and increase the over-all accuracy of the 
results by making the sweeping operations more consistent 
with the rest of the process. (3) Although the true wth 
eigenvalue is determined directly in the iterative transforma- 
tion procedure, the true nth mode must be computed from 
quantities within the iteration cycle after the transformed 
wth mode is found. 


OUTLINE OF STEPS IN THE PROCEDURE 

The equation of equilibrium of a cantilever beam vibrating 
harmonically in pure flexure is 


d?_ 

dx 2 


El 


d 2 y 
dx 2 


=yw 2 y 


( 1 ) 


or, after integration with proper attention to boundary 
conditions, 

( 'x ('x I f L n L 

V= Jo X ElJ, X (2) 

The solutions of this integral equation are the true natural 
modes (eigenfunctions) y l} y 2 , y 3 , . . . and the corresponding 
natural frequencies (eigenvalues) co I; w 2 , co 3 , . . . . For 
convenience in subsequent discussion, the true modes are 
assumed to be normalized to unity at some position (station 
A) along the beam. 

The first mode and frequency are assumed to have been 
previously determined by the Stodola process. The iterative 
transformation procedure becomes applicable in the deter- 
mination of the second mode and frequency. As mentioned 
previously, the immediate aim in the iteration for the second 
mode is the determination of a linear combination of first 
and second modes which is called the transformed second 
mode. The linear combination y 2 —yi which has zero 
ordinate at station A is chosen and defined as the transformed 
second mode to be determined. The iteration for determina- 
tion of this transformed second mode may be described as 
follows : 

(1) A plausible shape y u2 (1) for the transformed second mode 
is assumed. This shape must have zero ordinate at station 
A and should satisfy the boundary conditions as closely as 
possible. 

(2) The displacement 

J 'x pi 1 p L p L 

» X m XX 

resulting from the inertia load 7 co 2 2 y a 2 (I> corresponding to the 
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assumed shape y a 2 U) vibrating harmonically at frequency co 2 
is calculated. This calculation must usually be done 
numerically with the square of the frequency w 2 2 being 
carried along as an undetermined factor. 

(3) A first-mode shape (previously determined) is sub- 
tracted (swept out) from the calculated displacement y b in 
an amount such that the resulting displacement is zero at 
station A. Thus the resulting displacement is 

(4) The resulting displacement y a2 {2) is compared with the 
assumed displacement y a 2 (1) . When the computations are 
numerical, the ratios y a ^ 2) jy a 2 (1) are compared at all the sta- 
tions. If the assumed displacement is exactly equal to the 
transformed second mode, the ratios are equal to each other. 
These ratios contain the single unknown o> 2 , and the second 
frequency is that value of w 2 which makes the ratios unity. 

(5) If the ratios y a 2 l2) /y a 2 a) from the first cycle of iteration 
outlined in the four preceding steps are not reasonably the 
same at all stations, the process must be repeated until the 
ratios become reasonably the same. Each new cycle starts 
with the resultant displacement of each preceding cycle. 
The convergence of this process to the second frequency and 
the transformed second mode is proved in appendix A. 

The transformed third mode and the third frequency are 
computed in the following manner. The transformed third 
mode is defined as that combination of the first three natural 
modes which has a zero ordinate at the same station that was 
used in the transformed second mode (station A) and also a 
zero ordinate at some other station, station B. Thus the 
transformed third mode is defined as 

The iteration is as follows: 

(1) A plausible shape y a 3 (1) for the transformed third mode 
is assumed. This shape must have zero ordinates at sta- 
tions A and B and should satisfy the boundary conditions as 
closely as possible. 

(2) The displacement 


(4) The second sweeping operation, in which a transformed- 
second-mode shape (previously determined) is swept from 
the resulting displacement of step (3) so that the new 
resulting displacement is zero at station B as well as at 
station A, is performed. (This second sweeping operation 
cannot disturb the zero condition at station A established in 
step (3) because the second sweeping function (the trans- 
formed second mode) is identically zero at station A.) Thus, 
the final resulting displacement is 


(5) Comparisons of the ratios y a 3 (2 72/a3 (1) at all stations are 
made, and, if they are not reasonably the same, additional 
cycles of iteration are carried out until the ratios become 
reasonably the same. The third frequency is then computed 
from the ratios as explained previously. Convergence of 
this process to the third frequency and the transformed third 
mode is proved in appendix A. 

Frequencies and transformed modes higher than the third 
may be computed by extensions of the process just described. 

PHYSICAL INTERPRETATION OF THE PROCEDURE 

A physical interpretation of the iterative transformation 
procedure can be given. With regard to the transformed 
second mode, for example, the following question may be 
asked: Under what conditions can the beam vibrate in the 
transformed-second-mode shape at the second natural 
frequency? Vibration in shape y a2 — y 2 ~y l at frequency w 2 
implies an inertia loading yw 2 2 (y 2 — y\ )■ But if this load is 
substituted in place of yu 2 y in the right-hand side of equation 
(2), the result after integration will not be y 2 —y\ but 


V2 ~% V ^ £ L £ z w(?/2 -y^ dx v (3 > 

However, if an external (forcing) load of an amount 
7 (u 2 2 — w! 2 )y x is added to the inertia load, the total load 
y(u 2 2 yi— u 2 y\) will produce the shape y 2 —yi . Thus 



t(co 2 2 2 / 2 — <a 1 2 y 1 )(dx) i = y 2 — y i 


(4) 


»*= J. i El X X 

is calculated with the square of the frequency « 3 2 carried 
along as an undetermined factor. 

(3) The first of two sweeping operations, in which a first- 
mode shape is swept from the displacement y b so as to make 
the resulting displacement at station A zero, is performed. 
This operation gives the displacement 


2 / 6 — 



2/i 


The inertia and forcing loads are illustrated in figure 1. 
The inertia load acting alone produces a displacement 
(equation (3)) generally different from zero at station A. 
The forcing load produces the displacement 


(S - o w 


This displacement (equal to the sweeping function) has the 
shape of the previously determined first mode and is equal 
and opposite at station A to the displacement due to the 
inertia load; that is, by virtue of the previously assigned 
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normalizations at station A, 

Thus the displacement due to the forcing load is completely 
determined when the displacement due to the inertia load is 
known. The gist of the foregoing analysis is that vibration 
in the transformed-second-mode shape is the response of the 
beam to an oscillatory forcing load of the first-mode shape 
and of frequency equal to the second natural frequency, 
superimposed on a free vibration of the beam in the second 
natural mode. 

Similar physical interpretations of the iterative transforma- 
tion process for modes higher than the second can be made. 

APPLICATION OF THE PROCEDURE IN ORDINARY COUPLED 
NATURAL- VIBRATION PROBLEMS 

The procedure that has been outlined in a preceding section 
for pure flexure can easily be extended to systems capable of 
simultaneous flexural and torsional displacements. Airplane 
wings belong to the latter class of systems. The only 
distinguishing element in coupled flexural-torsional vibration 
problems is that each natural mode contains two components, 
the flexure and the torsion. These components must always 
appear together in a fixed relation to each other. The two 
components must be computed together and must be used 
together. 



(a) Transformed second mode: yai=yi—yu 

(b) Inertia load: TW2 s (!f2—ffi). 

(c) Forcing load: y(^2 2 —o>i 2 )yu 


Figure 1,— Illustration of physical basis of iterative transformation procedure. 


Each coupled mode is a solution of the simultaneous 


differential equations 

( 7 ) 

~^aJ^=W(Wy+ b V4!) ( 8 ) 

Equations (7) and (8), after integration, become (for a 
cantilever beam) 

?y== XX^X X ( 0 ) 

*^X (io) 


The solution of the integral equations (9) and (10) for the 
coupled transformed second mode by the iterative trans- 
formation procedure is outlined diagrammatically in figure 2. 
The flexural component of the displacement for a particular 
step is illustrated in the left-hand side of the figure and the 
torsional component is illustrated at the same level in the 
right-hand side. 



(a) Assumed transformed second mode. 

(b) Intermediate derived mode. 

(c) First-mode sweeping function. 

(d) Derived transformed second mode. 

Figure 2 . — Illustration of steps in the iterative transformation procedure for determining 

coupled modes. 
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In the first step, an approximation to a linear combination 
of the true first and second coupled modes is assumed. The 
particular linear combination having zero flexural displace- 
ment at the tip station (station A) is chosen. (For greatest 
numerical accuracy, this nodal point should be chosen in the 
component and at the station where the first coupled mode 
has its maximum numerical value.) The symbols y a2 (1) and 
<frx 2 (1) are used to designate the flexural and torsional com- 
ponents of this assumed displacement, respectively. In 
general, the magnitude of the torsional component relative to 
the flexural component is difficult to estimate ; the most exped- 
ient thing to do is to take one of the components equal to zero. 

The second step is the computation (by numerical inte- 
gration) of the two components of the displacement due 
to the inertia forces 7 « 2 2 (Va 2 -\-bu<j> a2 ) and inertia torques 
7 o> 2 2 ( & wy a 2 + 6 2 r 2 4 >a 2 ) that are associated with the assumed 
displacement. The result is termed the intermediate derived 
mode, and the symbols y b (l) and <£ 6 (1) are used to designate 
its two components. 

The third step is the determination of a sweeping function 
having the shape of the first coupled mode (previously de- 
termined) and a magnitude such that the sum of the inter- 
mediate derived mode and the sweeping function equals zero 
in the flexural component at station A. In algebraic terms, 
tbe first-mode sweeping function is given by 



(ii) 

---oa* 

(12) 


The fourth step is the addition of the intermediate derived 
mode and the first-mode sweeping function to give the derived 
transformed second mode. Thus the two components of the 
derived transformed second mode are 


yaa m =Vb w — 1 


(13) 


(£)/■ 

(14) 


The calculation of the ratios y a 2 ™/ya 2 W and ^ a2 <2) /0a2 (1> at 
all stations completes the first cycle of iteration. 

Additional cycles are carried out until the ratios at all 
stations in both the flexural and torsional components have 
values that are reasonably the same. The true second natural 
frequency of the coupled system is then obtained as described 
previously. 

Steady vibration of an airplane wing at zero airspeed is an 
example of coupled natural vibration. The actual numerical 
calculations for the transformed second mode as well as for 


the first mode and transformed third mode of an airplane 
wing vibrating at zero airspeed are discussed subsequently as 
a special case of flutter. 

The more general equations of airplane flutter at nonzero 
airspeed may be interpreted in such a way that they can 
be solved by a process analogous to that just described for 
coupled natural vibration. 

APPLICATION OF THE ITERATIVE TRANSFORMATION 
METHOD TO FLUTTER 

FLUTTER EQUATIONS 

The differential equations of equilibrium for a wing execut- 


ing simple harmonic motion are 

^ s Em+ig,)j^=-P,v+P t * (15) 

-^GJ(l+ig t )^-Q,V+Qt* ( 16 ) 

These equations govern a motion represented by 

Y(x,t)=y(x)e ial (17) 

$(x,t) = 4>{x)e iat (18) 


The use of the structural-damping coefficients g u and g$ in 
equations (15) and (16) is discussed in appendix B. The ex- 
pressions PyyArP^ and Q v y-\-Q^,<f> are the intensities of 
applied force and torque, respectively. For aerodynamic 
and inertia forces and torques due to air flow and distributed 
mass, the P and Q coefficients have values given by the follow- 


ing formulas (rearranged from those in reference 4) : 
For P y , 

Py — P Ry iPly 

in which 

(19) 

(20) 

and 

p -=(rJ 2 -fid^ 

(21) 

For P$, 

P<t>— Pr4> — iP til 

(22) 

in which 



Pr :* 


(23) 

and 


(24) 
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For Q v , 


Q,v — Qrv 1 Qiv 

(25) 

in which 


Qa= (lJ [-(!+“) x- a+tra ] (26) 

and 


Qn= (i) [~G +a )^ (£)>"’ 

(27) 

And for Qt, 


Q<t>— Qr <#> i Qi<t> 

(28) 

in which 


«*♦“(£) [ - (i _ ° i )u + G + ' i )f + § +<,!+ '‘’' 2 

E» v 

(29) 

and 


«»K0I(HHMS-G-)¥] 

( X ) b 0 W 
VM/o 

(30) 

For inertia forces and torques due to concentrated mass, the 
intensities of force and torque are, respectively, 

P,v+P t *= lim 

dx—*0 MX 

(31) 

and 

Q,y+Q t *= limfe'+erf 

dx — *0 ” X 

(32) 

in which 



(33) 


(34) 

and 



(35) 


For a cantilever wing the boundary conditions on the 
displacements are 

W,.o=(^.„=(g) i=o =[s/( 1 + ig.) 

=[^ 1+ ^SL A GJ(l+ig ^l-r° (36) 

The differential equations (15) and (16) are now written 
with the eigenvalue co 2 as an explicit factor. Thus equations 
(15) and (16) become 


and 


EI( 1 +ig.) g=^(P/»+P,» 

GJ{1 +ig,) =" ! «Vy+ QU) 


(37) 

(38) 


in which the P' and Q' coefficients are equal, respectively, to 
the P and Q coefficients divided by w 2 . 

FORMULATION OF PSEUDOFLUTTER PROBLEM 

Those solutions co 2 , ( y , <f > ) of equations (37) and (38) for 
which co 2 is a real and positive (not complex) quantity 
represent the steady harmonic motions of true flutter. 
However, because the P and Q coefficients are in general 
complex and because of the presence of structural damping, 
the solutions of equations (37) and (38) will, in general, 
be complex and will include complex eigenvalues co 2 . As 
in other methods of flutter analysis, the problem is made 
tractable by assuming at the beginning a value of the param- 
eter A:— — • This assumption fixes the values of the P and Q 
v 

coefficients. A real value of k is assumed because v must 
be real and only real values of co can represent flutter. To 
avoid the inconsistence of assumed real values of k and 
obtained complex values of co 2 in the solutions, the problem 
itself is altered by introducing an artificial damping so that 

co 2 

the complex eigenvalue is given by yq— r— ; where g a is the 

coefficient of artificial damping. Thus the differential 
equations of what may be termed the pseudoflutter problem 
become 

0 w «+«§-i 4 ( ww (39) 

WO) 

The value of co 2 can now be real for any assumed real value 
of k and is therefore the square of the frequency of the steady 
harmonic motion maintained by the artificial-damping 
forces and the naturally present aerodynamic, inertia, 
structural, and structural-damping forces. True flutter 
is possible for those special cases in which g a is zero. 

Equations (39) and (40) are similar in form to equations 
(7) and (8) and can be solved by the iterative transformation 
procedure in a way completely analogous to the solution 
of the ordinary problem. The complications introduced 
by the presence of air forces require, however, that a set of 
solutions be obtained for each of several assumed values of k. 
The fact that most of the functions involved are complex 
virtually quadruples the labor as compared with that 
required in the ordinary coupled natural-vibration problem. 

STEPS IN THE ITERATION AS APPLIED TO FLUTTER 

The iteration procedure employs the basic differential 
equations (39) and (40) in their integral forms which, for the 
cantilever wing under consideration, are 

(4i) 

*=M» z GJi^)Sp Q - y+Q/ * )(dxy (42) 


in which C stands for the more convenient form 


1 


of the 


eigenvalue. The iteration of equations (41) and (42) 
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follows the same form as the iteration of equations (9) and 
(10) . Briefly, the steps are as follows : 

(1) A real value of Jc is assumed and the values of the 
complex P and Q coefficients are computed. 

(2) An assumption is made for the desired mode y,<f>. 
(In the first cycle of iteration the assumed mode may be 
real but in the following cycles it will be complex.) 

(3) The complex loadings P v y-\-P$$ and Q, v y+ are 
computed. 

(4) The integrations indicated in the equations are 
carried out numerically to get the complex intermediate 
derived mode. 

(5) The sweeping operations are performed by using the 
complex lower-order transformed modes previously deter- 
mined. For convenience in numerical calculations, the 
flexural and torsional components of the complex derived 
(swept) transformed mode are computed in the forms 


and 


1 7o A 4 

C (Iq EqIq 


b 0 K v 


_1 To A 4 E 0 I 0 b q" tt- 
C no E 0 I 0 Go Jo A 2 ^ 


(43) 

(44) 


respectively, in which K v and K$ are nondimensional com- 
plex functions of the spanwise coordinate x. 

(6) The derived and assumed modes are compared by 
computing their ratios at the stations of the wing. If these 
ratios are not reasonably the same, additional cycles of 
iteration are earned out until the ratios are reasonably 
the same. In the limit (never obtained in practice) the 
ratios will be identical and the proper value of C is that value 
which makes them unity; that is, 

1 To A 4 , jr 1 To A 4 Eolo bo 2 jy. 

C Mo EqIq C mo EqIq GqJq L 2 ^ 1 

y 4> 


in which y and 4> constitute the assumed mode of the limiting 
cycle and the functions in the numerators constitute the 
derived mode of the limiting cycle. 

Equation (45) may be stated in the form 




To A 4 
Mo EqIq 


(46) 


in which D and H are nondimensional real numbers. Inas- 

1 % Q 

— ^2 ) equation (46) may be written 


much as C is defined as 


(D+iH) 


To A 4 ^ 1 +ig m 

Mo Eolo co 2 


(47) 


from which the frequency and artificial-damping coefficient 
are obtained as follows: 


V EqIq Hq 

■V 


9a jj 

217403 — 53 2 


(48) 

(49) 


The relative air velocity corresponding to the assumed value 
of k is given through the definition of k, that is, 



(50) 


NUMERICAL EXAMPLES 

Numerical computations presented in this section illustrate 
the actual application of the iterative transformation pro- 
cedure first to the ordinary natural-vibration problem 
(vibration at zero airspeed) and then to the flutter problem. 
All examples deal with the cantilever wing shown in figure 3. 

The geometric, structural, and mass properties of the 
wing are given in figure 3. A station coordinate system is 
employed for the purposes of the required numerical inte- 
grations. Four stations along the span have been selected 
as indicated in the figure; one of these stations is located at 
the spanwise position of the concentrated mass. The dis- 
tributions of forces and displacements over the span are 
considered to be adequately defined (through interpolation) 
by the forces and displacements at the four selected stations. 
The selection of a system of stations in any problem is im- 
portant because it greatly influences the amount and accu- 
racy of the work to follow. In problems, such as the present 
one, that involve concentrated masses, a station must be 
placed at each concentrated mass because displacements at 
the concentrated masses must be known. (More generally, 
a station must be placed at each discontinuity. Discontinu- 
ities may be present in the distribution of the structural 
stiffness and in the plan form as well as in distribution of 
mass.) The other stations should be equally spaced between 
the discontinuities, and for the system of parabolic inter- 
polation used in the numerical integrations in this report 
there must be a minimum of one station between each 
adjacent pair of discontinuities. The total number of 
stations should be the smallest possible that is consistent 
with the desired accuracy because the calculation effort 
increases rapidly with an increasing number of stations. In 
coupled systems, the number of degrees of freedom allowed 
is twice the number of stations selected; that is the number 
of degrees of freedom in either the flexural or torsional 



0.002378 slug/ft* 
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Figure 3. — Properties of cantilever wing used in numerical examples. 
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component of displacement is equal to the number of stations 
employed. Experience has indicated that with parabolic 
approximations results accurate to at least two significant 
figures in the highest mode computed can be obtained by 
employing numbers of stations as follows: For uncoupled 
systems, the number of stations should be two greater than 
the order of the highest mode to be computed; for coupled 
systems, the number of stations should be one greater than 
the order of the highest mode to be computed, with a mini- 
mum of three stations. More than these minimum numbers 
of stations may be required if their use is dictated by suffi- 
ciently many discontinuities. 

ORDINARY COUPLED NATURAL MODES AND FREQUENCIES 

The calculations for the first, second, and third modes at 
zero airspeed for the wing of figure 3 are shown in tables 1, 
2, and 3, respectively. In this case k= <» and the only 
aerodynamic forces are the apparent-mass forces. For 
simplicity, structural damping is disregarded; therefore, all 
quantities entering the problem are real. The numerical 
values of the aerodynamic-inertia force coefficients for 
k— co } as well as for other values of k to be used subsequently, 
are given in table 4. 

The first coupled mode is computed in table 1. Table 1 
shows in separate tabulations the flexural and torsional parts 
of the calculation. The first cycle of iteration (part (a) of 
the table) is shown in full detail. Two forms for the tor- 
sional part of the calculations are shown: The first form 

may be used when the torsional stiffness GJ is constant over 
each bay or over the whole length of the wing; the second 
form, which requires slightly more work, must be used 
when GJ is variable and may be used, as in this case, when 
GJ is constant. The second and third cycles of iteration 
are summarized in parts (b) and (c) of table 1. 

Details of the first cycle of iteration, if the procedure that 
applies only for constant torsional stiffness GJ for the tor- 
sional part of the calculation is used, are as follows: In 
columns 1 of table 1 (a) the two parts y x U) and <j> i (1) of the 
assumed first mode are listed. The torsional component is 
assumed to be zero because it will ultimately be small and is 
difficult to estimate. Columns 2 and 3 are the appropriate 
products of the assumed mode and the distributed-force 
coefficients, Columns 4, which are the sums of columns 2 
and 3, give the two components of the external load which 
correspond to the assumed mode and the arbitrary frequency w. 
Columns 5 give the concentrated loads (external forces and 
torques) that are equivalent to the distributed loads of 
columns 4. These equivalent concentrations are given in 
columns 5 in terms of the pertinent distances between sta- 
tions \i and in columns 6 in terms of the reference distance 
X 0 . Formulas used for computing the equivalent concen- 
trations from the distributed loads are given in appendix C. 
Columns 7 and 8 are the appropriate products of the assumed 
mode and the concentrated-force coefficients. Columns 9 
are the total concentrated loads, the sums of columns 6, 7, 
and 8. 


The flexural and torsional calculations must now be de- 
scribed separately. In column 10 for flexure, the average 
shears in the bays between stations are found by a successive 
summation of the concentrated loads from the tip where 
the shear is zero inboard to the root. In column 1 1 the incre- 
ments of bending moment are computed by multiplying 
the shears by the bay lengths in terms of A 0 . The bending 
moments of column 12 are found by a successive summation 
of the increments of bending moment from the tip where the 
bending moment is zero inboard to the root. Column 13 
gives the distribution of curvature, which is obtained by 
dividing each ordinate of the bending-moment curve by the 
local value of El (El in this example is constant) . Equiva- 
lent concentrated curvatures are now obtained by applying 
to the distributed curvatures the previously used formulas 
for equivalent concentrations. Column 14 gives these 
equivalent concentrations in terms of the distances A t -, and 
column 15 gives them in terms of the reference distance X 0 . 
The average slopes in the bays are obtained in column 16 
by a successive summation of the concentrated curvatures 
from the root where the slope is zero outboard to the tip. 
The increments of derived flexural displacement are com- 
puted in column 17 by multiplying the average slopes by the 
bay lengths in terms of X 0 . The flexural component t/i (2) of 
the derived mode is obtained in column 18 by a successive 
summation of the increments of displacement from the root 
where the displacement is zero outboard to the tip. Column 
19 gives the ratios at the selected stations of the derived 
flexural component to the assumed flexural component. 

Columns 10 to 15 for torsion are now considered. Column 
10 gives the average twisting moments in the bays of the 
wing and is obtained by a successive summation of the 
concentrated torques of column 9 from the tip where the 
twisting moment is zero inboard to the root. The average 
twists in the bays are computed in column 11 by dividing 
the average twisting moment in each bay by the local value 
of GJ (GJ in this example is constant over the whole span). 
The increments of derived torsional displacement are obtained 
in column 12 by multiplying the average twists by the bay 
lengths in terms of X 0 . The torsional component of the 
derived mode is computed in column 13 by a successive 
summation of the increments of displacement from the root 
where the displacement is zero outboard to the tip. Inas- 
much as the derived displacement of column 13 is in terms 
of GJ, the displacement is converted into terms of El in 
column 14 so that it may be compared with the assumed 
torsional displacement on the same basis as the assumed 
and derived flexural displacements are compared and so that 
the next cycle may be started with displacement components 
having the same dimensions as the assumed mode of this 
first cycle. Column 15 normally would contain the ratios 
at the selected stations of the derived torsional component 
to the assumed torsional component, but in the case of table 
1 (a) these ratios are meaningless because the torsional 
component will ultimately be different than was assumed 
in column 1. 
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Before the results of further cycles of iteration for the first 
mode are described, the form that the numerical integration 
for the torsional component must take when GJ is variable 
is described. In the part of table 1 (a) showing the calcula- 
tion for variable GJ, columns 1 to 4 are the same as in the 
calculation for constant GJ. The form of the numerical 
integration changes at column 5. Column 5 consists of 
increments of twisting moment over the bays. These incre- 
ments are obtained as increments of area beneath the curve 
of distributed torque (column 4). Formulas used for com- 
puting these increments are given in appendix C. In column 
5 the increments of twisting moment are given in terms 
of the distances Ai, and in column 6 they are given in terms 
of the reference distance A 0 . The twisting moments at the 
selected stations due to the distributed torsional loading are 
obtained in column 7 by a successive summation of the incre- 
ments of twisting moment. The components of external 
concentrated torque are as for constant GJ and are given in 
columns 8 and 9. The applied concentrated torque gives 
twisting moments as shown in column 10. Column 11 is 
the sum of columns 7 and 10 and gives the total twisting 
moments at the selected stations. (Note that in columns 10 
and 11 there is a discontinuity in twisting moment at the 
station having the mass discontinuity.) Column 12 gives 
the distribution of twist found by dividing column 1 1 by 
the local value of GJ {GJ being in general not constant). 
The increments of derived torsional displacement are com- 
puted in columns 13 and 14 by applying to the values of 
column 12 the same formulas applied previously to column 
4. The torsional component of the derived mode (columns 
15 and 16) is, except for small computational discrepancies, 
the same as in the previous method, as it should be. 

Two additional cycles of iteration were found to be ade- 
quate for the determination of the first mode and frequency. 
The results of these iterations are shown in parts (b) and (c) 
of table 1. In table 1 (b), for example, columns 1 give the 
two components of the assumed mode of the second cycle, 
which are obtained by normalizing the derived mode of the 
first cycle to unity in the flexural component at the tip 
station. This normalization is not essential but facilitates 
manipulations and comparisons by keeping the numerical 
values in all cycles within the same range of magnitude. 
Columns 2 give the derived mode obtained by the numerical 
integration procedure just described. The ratios of derived 
to assumed mode are given in columns 3 for both components 
of displacement. These ratios are seen to be fairly uniform. 
The ratios obtained in the third cycle in table 1(c) are, for 
practical purposes, identical. The averaging device shown 
in columns 4 of table 1 (c) and below table 1 is adopted 
as a quick and generally quite accurate way of smoothing 
out small discrepancies that remain in the ratios after 
convergence is almost complete. This device, although 
clearly not necessary in the case of table 1 (c), is useful in 
other cases throughout the numerical examples and is ex- 
plained as follows: The two ratios in columns 4 are obtained 
by considering the flexural and torsional components of the 


displacement separately and then dividing the sum of the 
station values of the derived displacement by the su m of the 
station values of the assumed displacement. When a dis- 
crepancy remains between two ratios of the type in columns 
4, the average of these two is taken as the final value; 
the final value for this case is given in the calculation 
below table 1 . This device gives greater weight to the larger 
ordinates and is in that respect similar to other weighting 
procedures such as the energy and least-squares methods 
but is much simpler. If the assumed and derived displace- 
ments contain both positive and negative ordinates, the 
negative ordinates should be changed to positive for the 
purpose of the summations. The remaining calculation 
shown below table 1 gives that value of the arbitrary fre- 
quency w which makes the ratio just computed unity. As 
proved in appendix A, this value of is the fundamental 
frequency «i. 

Table 2 gives the main results of three cycles of iteration 
required to obtain satisfactory approximations of the second 
frequency and the transformed second mode at zero airspeed 
{k— a a). Columns 1 of the first cycle (parts (a) of table 2) 
contain the two components y a2 (1) and <f> a2 (1) of the assumed 
transformed second mode. This mode must have one zero 
ordinate (excluding the root ordinates). Although this 
zero ordinate may theoretically be taken at any station, the 
numerical accuracy of the results is greatest if the zero ordi- 
nate is placed at the station and in the component where the 
preceding mode (the first) has its maximum numerical value 
(since the numerical process is such that the larger ordinates 
contain more significant figures than the smaller ordinates). 
Therefore, the zero ordinate of the transformed second mode 
is placed at the tip station in the flexural component, this 
location being designated station A. In the numerical 
values of columns 1, the flexural component y a 2 C1) would 
normally be taken as zero. (The values that are shown are 
estimated from flutter calculations that had previously been 
made for this wing.) Columns 2 give the intermediate 
derived mode obtained by numerical integration. Columns 
3 constitute the first-mode sweeping function. The shape 
of this sweeping function is given by columns 2 of table 1 (c) 
and its magnitude is such as to be equal and opposite to the 
intermediate derived mode at station A. Thus the derived 
transformed second mode (columns 4), which is the sum of 
columns 2 and 3, has zero ordinate in the flexural component 
at station A and a shape comparable to the assumed mode, as 
indicated by the ratios in columns 5. The ratios in the next 
two cycles (parts (b) and (c) of table 2) show marked im- 
provement in uniformity . The final value of the ratio com- 
puted below the table gives, as proved in appendix A, the 
value of the second frequency w 2 , as shown. 

The main results of the iterations to obtain satisfactory 
approximations of the third frequency and the transformed 
third mode at zero airspeed are stated in table 3. Typical 
operations required in a cycle are outlined in table 3 (a). 
Columns 1 give the assumed transformed third mode made 
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up of the two components y a3 l) and 0 a3 (1) . The transformed 
third mode is to have a zero ordinate in the flexural compo- 
nent at the tip station as in the transformed second mode and 
a zero ordinate in the torsional component at the tip station. 
The location of the second zero ordinate is designated station 
B. To obtain greatest numerical accuracy, the selection of 
the second zero ordinate is governed by the same rule that 
was used for selecting the first zero ordinate, namely, that 
the new zero ordinate should be placed at the station and 
in the component where the preceding mode (the transformed 
second) has its maximum numerical value. The numerical 
values that are shown in columns 1 are estimated from pre- 
vious flutter calculations; the torsional component <f> a3 w 
would normally be taken as zero. Columns 2 give the inter- 
mediate derived mode, and columns 3 give the first-mode 
sweeping function which, as before, has a magnitude at sta- 
tion A that is equal and opposite to the intermediate derived 
mode. Col um ns 4 constitute the transformed-second-mode 
sweeping function which has a shape given by columns 4 of 
table 2 (c) and a magnitude at station B equal and opposite 
to the sum of the intermediate derived mode and the first- 
mode sweeping function (the sum of columns 2 and 3). The 
derived transformed third mode of the first cycle is the sum 
of columns 2, 3, and 4 and is given in columns 5. The ratios 
in columns 6 are far from uniform. The ratios in the second 
and third cycles (parts (b) and (c) of table 3) show improve- 
ments in uniformity. The iteration is discontinued at the 
end of the third cycle where the ratios are about as uniform 
as they can get with the limited number of significant figures 
that are present. The frequency obtained by the smoothing 
device is the third frequency o> 3 and has the value shown. 

The patterns laid out in the foregoing examples establish 
the general technique that can be used to obtain zero- 
airspeed modes and frequencies higher than the third. 
Guiding rules for determining the number of selected stations 
to be employed have been given previously. These examples 
also set the basic pattern for the computation of the modes 
and eigenvalues of psuedoflutter and of flutter. 

MODES AND EIGENVALUES OF PSEUDOFLUTTER AND OF FLUTTER 

The operational solution in reference 5 gave for the wing 
under consideration (fig. 3) a reduced frequency at flutter 
of 0.1443. In order to use this operational solution, this 
same value (&=0.1443) is used in the flutter calculations 
that follow. 

The calculations for the first, second, and third modes at 
&= 0.1 443 are shown in tables 5, 6, and 7, respective^. 
Aerodynamic-inertia force coefficients have been computed 
by equations (19) to (35) and their values are given in 
table 4. Structural damping is disregarded, although a note 
on the method of incorporating structural damping in the 
calculations is made subsequently. 


Table 5 (a) shows in detail the first cycle of iteration for 
the first mode. The form of the computations is the same 
as that shown previously for the determination of zero- 
airspeed modes. The amount of computation, however, is 
between three and four times that required for zero-airspeed 
modes because of the fact that the functions involved are 
complex and thus must be described by two parts — a 
real part and an imaginary part. Columns 1 and 2 are the 
real and imaginary parts, respectively, of the assumed first 
mode. As a start, all parts of the assumed mode except the 
real part of the flexural component are taken as zero. Col- 
umns 3 to 6 are the real parts of the products of aerodynamic- 
inertia coefficients and the assumed mode, and thus their 
sums (columns 7) are the real parts of the distributed load. 
If the expressions for the distributed load are considered, 
this condition is more evident. The distributed forces pro- 
ducing flexure are given by 

(Pru — iP Iv) {y R + i Vi) + (P R<t> — iPl<t>) ( 4>R pi<Pl)=P Ryy R~\~ 

PRiti^RpPivyiP Pi4> < !>i J r'i(.PRvyi J rPR4>4 > i PjvVr Pi^r) 

(51) 

The terms of the real part of equation (51) appear in columns 
3 to 6 in the flexural part of table 5 (a); the terms of the 
imaginary part of equation (51) appear in columns 22 to 25 
in the flexural part of table 5 (a). This separation of real 
and imaginary parts allows the displacement due to each part 
to be computed separately. A similar explanation can be 
made for the quantities in columns 3 to 6 and columns 18 to 
21 in the torsional part of table 5 (a). 

Real and imaginary parts of the concentrated loads that 
are equivalent to the distributed loads are computed as 
explained previously by the formulas of appendix C. These 
values are shown in columns 8, 9, 27, and 28 in the flexural 
part and in columns 8, 9, 23, and 24 in the torsional part. 
The real and imaginary parts of the loads due to the con- 
centrated mass follow next in order, and the total concen- 
trated loads are given in columns 12 and 31 in the flexural 
part and in columns 12 and 27 in the torsional part. The 
average shears, average twisting moments, and bending 
moments are then computed as described previously. 

The remaining parts of the computations in table 5 (a) 
that are associated with the real parts of the load are de- 
scribed as follows (the remaining parts that are associated 
with the imaginary parts of the load are similar) : Column 
16 in the flexural part gives the distributed curvature due 
to the real part of the load. This curvature is obtained by 
dividing the ordinates of the real part of the bending-moment 
curve by the local values of the complex flexural stiffness 
El ( 1 + i <7 ?) ( 1 + ^<7 a ) ♦ In these examples, any actual struc- 
tural damping is disregarded; therefore g v is zero. The 
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factor l+i<7 a , containing the as yet unknown artificial- 
damping coefficient, combines with w 2 to give the factor 1/C 

in column 16, C being the arbitrary eigenvalue j— • 

If the actual structural damping g v is regarded as other than 
zero, the values in column 16 would be computed as follows: 
The real and imaginary parts of the bending moment would 
be combined into the complex bending moment 
This complex bending-moment distribution would then 
be divided by the local values of the complex stiffness to give 

R . ' wvv • — \* The factor 1 -\-iq a would be carried 
EI(l+ig y )(\+%g a ) 

along in the arbitrary eigenvalue C, and the numerical values 

of the real part of the quotient would be placed 

in column 16. The imaginary part of the quotient would 
be similarly placed (in column 35) in the calculations asso- 
ciated with the imaginary part of the load. The average 
twists due to the real part of the load are computed in 
column 14 in the torsional part of table 5 (a), and those 
due to the imaginary part of the load must also be computed. 
These calculations follow the same pattern as those just 
explained for the curvatures. The complex torsional stiff- 
ness GJ(l-\-ig^){l-\-ig 0 ) enters in place of the complex 
flexural stiffness. If GJ or is variable over a bay length 
or over the whole span, the numerical integration for the 
torsional part of the calculations should be carried out as 
explained in the part of table 1 (a) that deals with variable GJ. 

The numerical integrations are completed in the manner 
already described, and the derived mode is thereby obtained 
in the form of four components of displacement. The 
flexural components areyi B (2> and yu (2) of columns 21 and 40 
in the flexural part. The torsional components are <}> 1R l2) and 
<£u (2) of columns 17 and 32 in the torsional part. However, 
these components are not actually the real and imaginary 
parts of the flexural and torsional components of the derived 
mode, because each one of them contains the complex factor 
1+W7a- Nevertheless, the complex derived mode is given 
by yiR l2)J riyi/ 2) and -H<£u (2) . 

The complex ratios of the complex derived mode to the 
complex assumed mode are computed in column 41 in the 
flexural part and column 33 in the torsional part. Only 
two of these ratios have actually been computed but they 
are sufficient to indicate the need for further cycles of 
iteration. 


A total of four cycles of iteration (the main results of the 
last three are shown in parts (b), (c), and (d) of table 5) 
was required for satisfactory convergence. In columns 6 of 
table 5 (d) and immediately below table 5, the smoothing 
device described previously is applied to obtain the best 
single value of the ratios. The fundamental (first) eigen- 
value is that value of C which makes the ratio unity. Thus 


(7!= (269. 5 — 82.2-i) and since C x is defined as * Jr7 f al , 

JLL 1 jj. o)i 

the frequency and artificial damping of the first mode are 

obtained from the real and imaginary parts of the equation 

i+^l = (269.5-82.2i) (52) 

«i tLI fx 

The calculation of these quantities and the corresponding 
airspeed Vi which is obtained from the relation v i—~jr are 

shown at the bottom of table 5. 

Tables 6 and 7 show the main results of the iterations to 
obtain the transformed second and third modes for &=0.1443. 
Four cycles of iteration for each mode gave satisfactory 
convergence. The assumed modes of columns 1 and 2 of 
tables 6 (a) and 7 (a) were taken in the forms recommended 
previously in connection with tables 2 and 3. In tables 6 
and 7, the complex intermediate derived modes are given 
by ybR in) -\-iybi (n) and 4>b R in) +f</>w (n) , the complex first-mode 
sweeping functions, by y b \ R w -\-iy b u {n) and (j>bm w J ri<t>bu ln) 
with shapes corresponding to columns 3 and 4 of table 
5 (d), and the complex transformed-second-mode sweeping 
functions, by y ba2R w +iy b a2i w and <t> ba 2R W +i<l>b<t2i w with 
shapes corresponding to columns 7 and 8 of table 6 (d). The 
results computed in and below table 7 give for the third 
eigenvalue <7 a3 = 0.030 and co 3 — 168.9 radians per second. 
The corresponding airspeed is z? 3 =390 feet per second. 

COMPUTATION OF TRUE MODES 

Because the critical flutter velocities are given directly by 
the eigenvalues, knowledge of the true modes in flutter 
problems is of no value (at least of no value recognized at 
present). The same statement applies to the transformed 
pseudoflutter modes, with the exception that in the iterative 
method their determination is a necessary adjunct to the 
determination of the eigenvalues. In ordinary problems of 
forced vibration (at zero airspeed), however, the true modes 
are often used with great advantage. For this reason and 
for the sake of completeness of the presentation of the itera- 
tive transformation procedure, the method of determining 
true modes from results of the iterative transformation 
procedure is illustrated in tables 8 and 9. 

The computations in tables 8 and 9 pertain to the same 
wing analyzed in the previous examples. The modes com- 
puted are for &=0.1443. The true third mode as computed 
in table 9 may therefore be compared with the flutter 
mode computed for this wing by the operational method in 
reference 5. 

In table 8, the true second mode is computed as follows 
from functions appearing in the last cycle of iteration for the 
transformed second mode (table 6(d)): Preceding the table 
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c 

proper is the calculation of the eigenvalue factor F 12 =jj — 1 


that is needed for computing the true second mode. In the 
terminology of tables 6(d) and 8 and as shown in appendix A, 
the true-second-mode shape is given by 



y 2R + iya= (yiR+iyu ) + 4-^W 5) ) 

(53) 

and 

<t>2R + i <t>21 = (4>1R 4~ 4" (<f>02 R (5> + ^ <W ( ® ) 

(54) 

in which 

1 n*+m,- v '‘* W c +iy ™ W 

(55) 


cr 1 


and 

. , <f>61R (4) +^6U (4) 

1 Q 

(56) 





Columns 1 of table 8 show the key ordinate (t/ S i B c4) +iytu w )A of 
the first-mode sweeping function y blR (i) +m/mj (4) , 4>bm w +i<f>M/ (4> 
as given in columns 5 and 6 of table 6(d). The key 
ordinate is taken as the largest ordinate (the ordinate at 
station A) for the reason of accuracy cited previously. The 
key ordinate of the first-mode shape yw+iyu, 4>\ R+i4>u (equal 
to the first terms on the right-hand sides of equations (53) 
and (54)) is shown as the boxed value in column 2 and is 
obtained by dividing the value in column 1 by the eigenvalue 
factor F 12 . The other values in columns 2 are obtained by 
using the key ordinate in conjunction with the first-mode 
shape given in columns 3 and 4 of table 5(d). Columns 3 
show the transformed-second-mode shape y a 2 R (S> +«/a2/ <5) , 
(j>a 2 R (S) +^a2r (5) (equal to the second terms on the right- 
hand sides of equations (53) and (54)) given by columns 
7 and 8 of table 6(d). The sum of columns 2 and 3 which 
is given in columns 4 gives the shape of the true second 
mode y 2 R+iy 2 i, (equal to the left-hand sides of 

equations (53) and (54)). 

In table 9, the computation of the true third mode pro- 


a 


ceeds as follows: The necessary eigenvalue factors F l3 =-~-— 1 

'- / 3 

Ct 

c 3 


c 

and I^3=7T — 1 are computed as shown. In the terminology 


of tables 7(d) and 9 and as shown in appendix A, the truc- 
third-mode shape is given by 


y^R + i y 3 / = (y ur 4 - i ym) + (y i?r + i y 127) 4 - {yaw + iy a 2i) 4 - 


(W 5) +W 5) ) 

and 


( 57 ) 


4>zR-\-i4>zi~{^uR J ri<f}ui) J r{.4>\2R J r‘i4>m)- i r{^a2R J ri<i>a2i) J r 


(</>«3/J (5) +?>a3J (5) ) 


(58) 


in which 


(^JJ~ l^( 2 /niz+'iyn/)+^£T — - 1 ^) (y ur^t iy 121) 


=y bl R w 

+ iybu w 

(59) 

f ^ (<#> 11 ^ + ^ 117 ) — 1^) (<l>l2R~\~'i' < f>12l) 


-C> 

-e- 

II 

+ 'f<f>617 (4) 

(60) 

ya2R~\~ ^2/a27‘ 

y b a2R W +iyba2I W 

(61) 

F 3 



4 > a2R~\~ 1 4 > a2I : 

4>ba2R {i) Fi4>ba2I < ' i) 

(62) 

c 2 


c 3 1 



and y 12R +iy 12 r, <Pi 2 R-\-i<t>i 2 i is to y a2 R+iy a2 r, baiR+i^i as 
ViK+iyu, is to y a 2R (r,) +iy a 2i (5) , <t>a 2 R (5) -K<w 5) in table 

8. The key ordinates (y bl R W +^2/m/ (4) )a and (<j> ba2 R (i) J ri4>M2i {A) )B 
of the first and second sweeping functions appear in columns 1 
and 2 and are taken from columns 5, 6, 7, and 8 of table 
7(d). The key ordinate of the functions F 23 (y 12R -{-iyi 2 i), 
F 23 {4>\ 2 R-\-i<i>i 2 t)> which are equal to the second terms on the 
left-hand sides of equations (59) and (60), is computed in 
columns 3 by using the key ordinate of columns 2 in con- 
junction with ordinates at stations A and B in columns 2 
and 3 of table 8 as follows: 


[F 23 (y V2R + ^ 2 / 12 /) a] table 9 

{yiR-\-iyu)A 


> 


(<w (5 ) +^<W 5 ) )bJ table 


1 


\F 23 {4>a2R 4 " i </>a 2 /)s] table 9 ( 63 ) 


The key ordinate of the functions ^3(^1^ -f iy lu ), F l3 (<f> nR -\-i<i>iu), 
which are equal to the first terms on the left-hand sides 
of equations (59) and (60), is given in columns 4 and, 
in accordance with equations (59) and (60), is the difference 
between y blR (i) 44y&u (4) , </>&ib ( 4) +uf>&i/ (4) of columns 1 and 
F 23 (y 12B +iy 12 i), F 23 ((f>i 2 R-\~i<t>i 2 i) of columns 3. The key ordi- 
nates of the first-mode shapes yn/j + M/n/, and 

yi 2 R-\~iyi 2 i, </ , i2K+'i</>i27 are shown in columns 6 and 5 and are 
obtained by dividing the values in columns 3 and 4 by the 
appropriate eigenvalue factors. The sum of the key ordi- 
nates of columns 5 and 6, shown as the boxed value in col- 
umn 7, is the key ordinate of the total-first-mode shape 
2/ifl+M/m <t>iR-\-i<l>u which is equal to the sums of the first two 
terms on the right-hand sides of equations (57) and (58). The 
other values in columns 7 are obtained by using the key ordi- 
nate in conjunction with the first-mode shape given in columns 
3 and 4 of table 5 (d) . The key ordinate of the transformed- 
second-mode shape y a2 R J riya 2 i, <?W44<?W, which is equal 
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to the third terms on. the right-hand sides of equations 
(57) and (58), is shown as the boxed value in columns 8 and is 
obtained by dividing the value in columns 2 by the eigen- 
value factor F 2S . The other values in columns 8 are com- 
puted by using the key ordinate in conjunction with the 
transformed-second-mode shape given in columns 7 and 8 of 
table 6(d). Columns 9 show the transformed-third-mode 
shape y M m +iy a3 i w , <?W 5) +^<W s> (equal to the fourth 
terms on the right-hand sides of equations (57) and (58)) 
given by columns 9 and 10 of table 7(d). The sum of col- 
umns 7, 8, and 9 given in columns 10 gives the shape of 
the true third mode yzR-\-iy 3 i, (equal to the left- 

hand sides of equations (57) and (58)). 

TRENDS AND COMPARISONS OF NUMERICAL RESULTS 

Results of the computations shown in the preceding 
section of the report together with results of similar compu- 
tations based on other assumed values of k are given in 
figures 4 to 6. Figures 4 and 5 deal with the wing to which 
the concentrated mass is attached. Figure 6 gives data of a 
s imil ar nature for the same wing without the concentrated 
mass. The computed results obtained by the Rayleigh- 
Ritz and operational methods and the experimental results, 
all of which are given for this wing in references 5 and 6, are 
also recorded in figures 4 to 6. 

In part (a) of figure 4 the solid curves show the variation 
of the artificial-damping coefficient g a with airspeed in each 
of the first three solutions. For each assumed value of k a 
dashed curve is drawn through points that represent solutions 
for that value of k. Part (b) of figure 4 shows in a similar 
way the variation of the frequency a with airspeed and the 
lines of constant values of k. The facts of particular interest 
that are shown by these plots are as follows: 

(1) The true flutter condition is given by the third solution 
for a value of k between 0.1443 and 0.1590 at an airspeed 
almost equal to that found in the experiment. Here the 
computed value of g a is zero. The computed frequency at 
true flutter is also in very close agreement with the experi- 
mental value. 

(2) The operational solution is in good agreement with the 
experimental solution, but the solutions obtained by the 
Rayleigh-Ritz method v T ith three and four modes vary by 
72 percent and 22 percent, respectively, from the solution 
obtained by the operational method. The operational solu- 
tion is theoretically the most exact even though it involves 
summations of finite numbers of terms of infinite series. 
However, as pointed out in reference 5, its use is limited in 
practice to wings of uniform section. In the present example 
the results obtained by the iterative method would be 
expected to be better than the results obtained by the 


Rayleigh-Ritz method because the eight degrees of freedom 
used in the iterative method are much less restrictive than 
the three or four used in the Rayleigh-Ritz method. 
Although exact agreement of the results of any of the com- 
putational methods with the experimental results is not to 
be expected, the better agreement of the iterative solution 
as compared with the operational solution is at first sur- 
prising. On further observation, however, this agreement 
must be credited to a fortunate disposition of the errors 
involved in the -iterative method because, in the case of 
figure 6, the relative order of agreement of the operational 
and iterative results with the experimental result is opposite 
to that in figure 4. 

(3) The trends of the solid curves representing the first 
and second solutions in figures 4 (a) indicate that both may 
cross the zero artificial-damping axis at very large airspeeds. 

0 Experiment (reference 5) 

Operational solution (reference 5) 

* Rayieigh-Ritz •• 3 modes (reference 6) 
v Rayieigh-Ritz 4 modes (reference 6) 

° iteration: 4 stations 



(a) Variation of artificial damping with airspeed. 

(b) Variation of frequency with airspeed. 

Figure 4.— Variation of artificial damping and frequency with airspeed in first three solutions. 
Wing with concentrated mass, 
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But this conjecture is of no practical interest so long as a 
curve (the third solution) that crosses at a lower airspeed 
exists. However, the question of whether the curve for 
some solution higher than the third could cross the zero 
artificial-damping axis at an airspeed lower than that at 
which the third solution crosses demands an answer. 

(4) Reasonable assurance that, among all possible solu- 
tions, the curve of third solutions in figure 4 (a) crosses the 
zero artificial-damping axis at a lower airspeed than any 
other is provided by the trends of the curves for constant 
values of £ in parts (a) and (b) of figure 4. The curves of £ 
show that the curve representing the fourth solution will 
most assuredly lie above and to the right of the solid curves 
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(a) Variation of artificial damping with airspeed. 

(b) Variation of frequency with airspeed. 

Figure 6. — Variation of artificial damping and frequency in first three solutions. W ing 
without concentrated mass. 


in figure 4 (b) and probably below and to the right of the 
solid curve for the third solution in figure 4 (a) . The curves 


of it in figure 4 (b) are straight lines by definition 



Prediction of the courses of the curves of £ in figure 4 (a) 
cannot be made with much certainty. They have a strong 
tendency to proceed to the right, but it is easy to believe that 
upward or downward changes in their directions could take 
place. The curve for the fourth solution, however, would 
probably cross the zero damping axis at a value of v between 
500 and 600 feet per second in figure 4 (a) . 

Figure 5 shows and compares the amplitude and phase 
distributions of modes computed by the iterative transfor- 
mation procedure and by the operational method for the 
wing with a concentrated mass. The first and second modes 
as well as the more important third mode from the iterative 
solution for k= 0.1443 are plotted, and the third mode from 
the iterative solution for £=0.1590 is also plotted. The 
third modes from the iterative solutions for the two values 
of k agree very well in shape with the flutter mode obtained 
in reference 5 by the operational method, and the operational 
mode lies between the two iterative modes. Thus the 
agreement of the iterative and the operational methods is 
again evidenced. 

Figure 6 is a plot similar to figure 4 but relates to the 
behavior of the wing analyzed in figure 4 if the concentrated 
mass is not present. There is very little similarity in the 
data of the two figures. The most notable difference is that 
in figure 6 the true flutter mode appears in the second solution 
instead of the third as before and that the flutter speed is 
lower than before. Of interest is the occurrence of almost 
equal eigenvalues in the second and third solutions for 
£=0.50. The flutter speeds given in figure 6 by all methods 
of solution, including the Rayleigh-Ritz method, are seen to 
be in substantial agreement. 


CONCLUDING REMARKS 


The report has described the iterative transformation 
method suggested by H. Wielandt and has demonstrated 
the use of the method in an orderly computation of critical 
flutter speeds. Numerical comparisons with solutions ob- 
tained by other methods and with experimental values have 
been made. The applications made in this report show 
promise for future practical use of the method. 


Langley Aeronautical Laboratory, 

National Advisory Committee for Aeronautics, 
Langley Field, Va., January 17, 1951. 
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APPENDIX A 

ON THE CONVERGENCE OF THE ITERATIVE TRANSFORMATION PROCEDURE 


INTRODUCTION 

The extensive existing literature on the eigenvalue prob- 
lems is concerned almost exclusively with the class known as 
self-adjoint problems, in which the eigenfunctions and eigen- 
values are real. In recent years, non-self -adjoint eigenvalue 
problems have received increasing attention. This class 
includes the flutter problem in which the eigenfunctions and 
eigenvalues are generally complex. The literature referred 
to by Wielandt in reference 3 reveals that the non -self-adjoint 
eigenvalue problem and the transformation method for its 
solution have been given some attention since at least 1928. 
Wielandt’s own work constitutes probably the most extensive 
contribution on the subject. 

The discussion on convergence given herein is not con- 
tained in Wielandt’s work and ma} r be considered a rigorous 
proof if the following assumption is valid: that the equations 
(equations (41) and (42)) for the system (the wing) under 
consideration have an infinite number of solutions that form 
a complete set for any value of the reduced frequency Jc. 
In the subsequent demonstrations, the validity of expanding 
arbitrary displacement functions in infinite series of eigen- 
functions depends upon the validity of the assumption. 
That complete sets of eigenfunctions do exist seems plausible 
enough to justify reliance in the conclusions. 

BASIC RELATIONS 

For any one of the true solutions of the eigenvalue problem, 
for example, the eigenvalue C m and eigenfunction y m ,4>m, 
equations (41) and (42) may be written as 

Jo EI(l+ig v ) X X (P v 'y m ~hP/<P m )(dxy (Al) 

and 

= Gj{\+"igt) § x (Qv'ym + Q*'<i>m)(dxy (A2) 

To make the notation more concise, let the coupled mode 
y m ,<f> m be represented by w m . Then if y m ,<i>m. is substituted 
into the right-hand sides of equations (Al) and (A2), the 
left-hand sides may be represented by C m w m ■ Furthermore, 
because of the linear character of the equations of the prob- 
lem, substitution of the function series 

(A3) 

i = 1 

into the right-hand sides of equations (Al) and (A2) gives 
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for the left-hand sides the function series 

X ) CiCLiWi (A4) 

i=l 

The coefficients a { are, in general, complex. The complex 
eigenvalues Ci are assumed in the subsequent proofs, except 
where stated otherwise, to be different from each other, and 
the eigenvalue having the largest modulus is defined as C\, 
the second largest, as C 2 , and so forth, so that 

ic,i>|c,|>icy> . . . (As) 

Expressions (A3) and (A4) are the expansions, in terms of 
the eigenfunctions and the eigenvalues, of the functions 
previously referred to as the assumed and intermediate de- 
rived modes, respectively. The subsequent proofs of con- 
vergence are based upon the fundamental relationship that 
exists between expressions (A3) and (A4). 

FUNDAMENTAL MODE 

The fundamental mode and eigenvalue are found by itera- 
tion according to the original Stodola procedure. In the 
present terminology and notation, this procedure and its 
proof of convergence are as follows: The coupled mode as- 
sumed at the beginning of the first cycle of iteration in general 
contains some component of each of the eigenfunctions; 
therefore its most general expression is 

Wi (1) =X) a t w i (A6) 

The intermediate derived mode (which in this case is also the 
final derived mode inasmuch as no sweeping operation is 
required to obtain the first mode) is for this first cycle of 
iteration 

w b (1) = w 1 (2) =2 OidiWi (A7) 

; = i 

The second and following cycles are begun with the final 
derived mode of each preceding cycle, and thus the assumed 
and derived modes of the nth cycles are 


co 


1= 1 

(AS) 

Wi (n+1) =S Cfam 

i— 1 

(A9) 
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In accordance with the definitions given in equation (A5), 
all terms on the right-hand sides of equations (A8) and (A9) 
except the first are negligibly small in comparison with the 
first for large values of n. In the limit the fundamental mode 
is obtained as 

lim w 1 (n+1) =lim Ci n aiW x (AlO) 

71— »co 71— » co 

and the fundamental eigenvalue is obtained from 

w (ra+l) 

lim ^k*r=°' (An) 

ft— » 0 D M/l 

TRANSFORMED SECOND MODE 


When each succeeding cycle is begun with the derived 
transformed second mode of its preceding cycle, the various 
functions for the 7ith cycle are 


Wa2 (B) = S C i n ~ l b i {Wi — W i) 
1=2 

(A17) 

W6 (n) = ib O i n ~ l b lip iWi—C jWj) 

i= 2 

(A 18) 

w fl2 (n+1) — S Ci n bi(Wi—Wi) 
1 = 2 

(A19) 


The limits as n approaches infinity are 


The initial assumption of the transformed second mode 
in general is of the form 

<A12) 


and 


lim w a2 (n+1) =lim C 2 n b 2 (w 2 — w x ) 


lim 

71 — » co 


w a2 (n+1 > 

w a2 w 


C 2 


(A 20) 
(A21) 


in which the arbitrary coefficients bi are in general complex 
and the subscript A refers to values of either the flexural or 
torsional components of the eigenfunctions at station A. 
More specifically, if, for example, the nodal (zero) point of 
w a2 is selected to be at station A in the flexural component, 
then the subscript A refers only to the flexural components 
of Wi, w 2 , w 3 , . . . and not to their torsional components. 
Thus each term of the series in equation (Al2) satisfies the 
requirement that either the flexural or torsional component 
of the assumed mode be zero at station A. 

To simplify the subsequent work as much as possible, the 
eigenfunctions are henceforth assumed to be normalized to 
unity at station A; thus 


Equations (A20) and (A21) show that convergence to the 
exact-transformed-second-mode shape vu 2 —w x and to the 
exact second eigenvalue C 2 can be obtained theoretically. 

TRUE SECOND MODE 

The key to computation of the true second mode is readily 
found in the simple case illustrated in figure 1. In this case 
the sweeping function of the final cycle of iteration would be 
the displacement produced by the forcing load y(w 2 2 — coji 2 )^, 
in which y x is the first-mode component of the transformed 
second mode y a2 . The sweeping function is designated by 
y si which has a well-defined numerical value in the iteration. 
Thus the value of y x could be found from the equation 


Wa= 1 (i= 1,2,3,...) (A 13) 

Equation (A12) now takes the simpler form 

w a2 (1) =ib bi(Wi—Wi) (A 14) 

i=2 

The assumed mode given by equation (Al4) leads, ac- 
cording to equations (A3) and (A4), to the following inter- 
mediate derived mode: 


that is, 


2/m = 


7(co 2 2 — -cji 2 )y x 
7cV 



2/1 


2/61 

o) 2 2 

«1 2 


(A22) 


(A23) 


The sum of y a2 given in the iteration and y x given by equation 
(A23) gives y 2 , the true second mode; that is, 


w 6 (1) = X) bi(GiWi—CiWi) (A 15) 

i=2 

Sweeping of this intermediate derived mode with the first- 
mode shape (previously determined) leads to the derived 
transformed second mode of the first cycle as follows: 

Wl = ±; C t b t (w t -w^ (A 16) 

\ J A i=2 


2/«2 + 2/i= 2/2— 2/i+2/i = 2/2 (A24) 

By analogy, the true-second-mode shape in the general 
(complex) problem under consideration is found as follows : 
The limiting value of the sweeping function is, from equation 
(A18), 

lim w 61 (n) =lim ir^lim C 2 n (yr— 1^) b 2 w x (A25) 

«— * co n—*<D \ W i /a. n— ><= \U’2 / 
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The expression analogous to equation (A23) is 
lim w b i (n) 


=lim C 2 n b 2 Wi 


Cl ■, n— 

c 2 

The expression analogous to equation (A24) is 

lim w 61 (n) 


lim Wa 2 (n+1) ■ 

'll — » CO 


=lim C 2 n b 2 w 2 

Oj 1 71— >co 

C 2 


(A 2 6) 


(A27) 


which gives the exact shape of the true second mode. 


TRANSFORMED THIRD MODE 


Sweeping of the mode given by equation (A31) with a 
transformed-second-mode shape such as to make the sum zero 
in the flexural or torsional component (as the case may be) at 
station B gives the derived transformed third mode as follows: 


W aZ 


/w b in) \ 

| Wi — 
A 

W b (n) — I 

\ Wi ) 

) Wy 

A | 

\ Wi ) 

L w 2 —w 1 

__ B 


1 =W 6 (n) 

= (w 2 —w x ) 1 (A32) 

i ~ 3 L \W 2 ~Wi/ B J 


The limits as n approaches infinity are 
lim w a3 (n+1) =lim C*d z fw 3 — w x — ^ (w 2 ~w{) 1 (A33) 

n — >“> n — L \W 2 W\J b _J 


The first cycle of iteration for the transformed third mode 
begins with an assumed mode that has two zero values, one 
of these being in the same (flexural or torsional) component 
and at the same station (station A) as previously employed 
for the transformed second mode. The other zero value 
may be taken in the same component as was the first zero 
value and at a different station (station B), or it may be 
taken in the other component at any station, including 
station A. Either of these possible selections for the loca- 
tion of the second zero value is indicated in the following 
equations by use of the subscript B. The initially assumed 
transformed third mode may be written as 

M>.. <u = Srf«r (w.-».)l (A28) 

i=3 L \Wi — Wi/ B J 

in which the arbitrary coefficients d* are complex. Each 
term of the series in equation (A28) is zero at station A by 
reason of the normalizations stated in equation (Al3), and 
each term is also zero at station B. 

The various displacement functions for the general (wth) 
cycle of iteration may be expressed as follows: The assumed 
mode is 


w^"<=± (»,-«,)! (A29) 

i= 3 L \ W 2 — Wi/b J 

The intermediate derived mode is 

w„^ = ± {C 2 w 2 -Ciw x ) 1 

i=3 L \w 2 — Wi/b J 

(A 30) 

The result after sweeping the intermediate derived mode 
with a first-mode shape such as to make the sum zero at 
station A is as follows: 


w 


C n ) 


Wl =£ C t *-'dS CM-wd- 

V wi Ja m L 

C 2 (w 2 - Wl ) 1 
\w 2 —wJ B 'J 


and 

vi 

lim "w ~w ~ =g 3 (A3 4) 

As shown by equations (A33) and (A34), convergence to the ex- 

act-transformed-third-mode shape w 3 — Wi — (w 2 —w 1 ) 

and to the exact third eigenvalue C z can be obtained theoreti- 
cally. 


TRUE THIRD MODE 

Computation of the true third mode is explained by refer- 
ring again to the simple problem of pure flexural vibration 
in which air forces are excluded. The transformed third 
mode in this simple problem would be given by 

yas=ys-yi~(~~) — yO (A35) 

\y 2 —yi/ B 

The total load required to hold the beam in equilibrium in the 
shape y aZ is 

[i J ( A36 > 


If the beam is vibrating with shape y a3 at frequency co 3 , the 
inertia load is given by 

7u 3 2 ^3=7w 3 2 Q/3— y\— (A37) 

The forcing load required is the difference between the total 
load (expression (A36)) and the inertia load (equation (A37)), 
that is, 




(A 3 8) 


The displacement produced by this forcing load is 


A(31) 
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and this displacement must be equal to the sum of the 
sweeping functions in the last cycle of iteration (if the 
iteration has been carried to complete convergence). The 
first sweeping function is of the first-mode shape and the 
second sweeping function is of the transformed-second-mode 
shape. If the expression (A39) is written in the form 


each of the sweeping functions contained in the displacement 
produced by the forcing load is obvious. Thus 




and 




(A41) 

(A42) 


in which y b i and y ba2 designate the first and second sweeping 
functions, respectively. Both of these functions have well- 
defined numerical values in the iteration. 

If now a simpler notation is adopted, equations (A41) and 
(A42) can be written as 




(A43) 

and 

(A 4 4) 

v “ 

in which 

(A45) 


and 


(A46) 



(A 4 7) 

The true third mode is clearly given by the sum 
(A35), (A45), (A46), and (A47); thus 

of equations 


Vz — y<ia 4 - yi 1 4- y 1 2 + yea. 

(A48) 


The transformed third mode y a 3 is given directly in the 
iteration. The procedure for finding the other components 
on the right-hand side of equation (A48) is as follows: 
Component y a2 , by equation (A44), is 

(A 49 ) 

9 -L 


Component y n is known when y a i is known because its rela- 
tion to y a2 was established previously in connection with the 
transformed-second-mode calculations (see equation (A24)). 
Component y n is then found by equation (A43) as 



By analogy with the foregoing case, the true third mode in 
the complex-eigenvalue problem is found as follows: The 
limiting value of the second sweeping function is (see equa- 
tions (A31) and (A32)) 


lim w ba 2 W = 


ro ,i 

— lim V k .. 

ft— *co — *10 X 1 


(w 2 —w 1 ) 


-lim C 3 n (%- 1 (w a — wO (A51) 

■n— >0) \ 0 3 / \®2 Wt/B 


The limiting value of the first sweeping function is (see 
equations (A30) and (A31)) 


lim w bl (n) = — lim (— \ w x 

n — >03 n — >« V Wi /a 


=lim C 3 n d 3 


Vi_ 1 ( w 3 — wA “1 

C 3 \C 9 Cj \w 2 — wJbJ 


Wi 


(A52) 


The quantities analogous to y i2 and y a 2 of equations (A43) 
and (A44) are, for the present case, lim Wi 2 (n) and lim w a2 in) . 

ft— >00 ft— >co 

The latter quantity is obtained from the relation analogous 
to equation (A49) as follows: 


lim w a2 in) 

n—> ® 


lim w ba 2 < ' n) 

ft— » CO 

c 3 


dim C 3 n d 3 C~—^\ (w 2 -w{) (A53) 

ft— >00 \ wi / B 


The relationship of lim Wi 2 (n) and lim w a2 w is obtained from 

ft— >00 ft— » co 

equations (A26) and (A20) of the section dealing with the 
transformed second mode. Thus, 


lim 

ft— > co 


w n (n) —lim 

ft— >co 


[{Wa 2 in) ) B ] Eg. C/153) 
[(W a2 (n+1) )B] Eg. (.420) 


lim w b i w 


ft— > CO 



Eq. (.426) 


=lim C 3 n d 3 



(A54) 


The quantity analogous to y n of equation (A43) is, for the 
present case, lim w n <n) and is obtained by an equation 

ft — >00 


20 


REPORT 1073 — NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS 


analogous to equation (A50) as follows: 

lim w bl w A lim Wi 2 (n> 

lim ^ 

n-*°° x .. 

c t 

=lim C 3 n d 3 f l—(— — — ^ 1 Wi (A55) 
n-*» L \w 2 — wjsj 

The exact shape of the true third mode w 3 is given by the 
sum of equations (A33), (A53), (A55), and (A54), which is 

lim {w a3 (n+1) +w o2 (n) +w n (n) + ^ 12 <re) ) =lim C 3 n d 3 w 3 (A56) 

n—> oo n— 

FOURTH AND HIGHER MODES 

Extensions of the proofs to modes higher than the third 
can be made in a manner similar to the foregoing proofs. By 
this means, the iterative transformation procedure can be 
proved, under the assumptions stated at the beginning of this 
appendix, to be convergent for all modes and eigenvalues. 

CASES OF EIGENVALUES HAVING EQUAL OR NEARLY EQUAL MODULI 

For a representative case, suppose that 


and derived modes in each of the cycles contain only com- 
ponents of the types in equations (A60) and (A61). 

The two components (with shapes w 2 —vji and w 3 —Wi) 
clearly appear in the last cycle in proportions different than 
in the preceding cycle. (The proportion in each cycle is a 
complex function of the spanwise coordinate.) Because of 
this differing proportionality the results of cycles n — 1 and n 
can be linearly combined so that the combined functions 
contain only one of the components w 2 — Wi and w 3 —Wi. 
Accordingly, the ratios of both the flexural and torsional 
components of the combined functions at all stations should 
be equal to each other. In algebraic terms, this statement 
means that 


( rw a2 {n) +w a2 ( " +1) \ 
\rw a 2 in ~ l) +Wa2 (n) )s 


(A62) 


in which r and R are (complex) constants, and the subscript 
S designates that the ratio may be evaluated at any station 
S, that is, that R has the same value for all stations. All w 
functions must be the same type of component, either flex- 
ural or torsional. 

Since S can be any station, the equality 



\Ci\>\C 2 \i l<? 3 |>|U 4 |> . . . 

(A57) 

and that 


icy = icy 

(A58) 

or that 


\C 2 \ ~ \C 3 \ 

(A59) 


Under conditions (A57) and either (A58) or (A59), the 
assumed and derived modes after a few cycles of iteration 
will be virtually as follows (see equations (Al7) and (Al9)): 

w a2 (,l) = C 2 n ~ l b 2 (w 2 — w0 + 0^-% (w 3 — Wj) (A 60) 

w a2 (n +1) = C 2 n b 2 ( w 2 — w?i) + C 3 n b 3 (w 3 ~ wO (A 6 1) 

If \C 2 \ is only slightly greater than \C 3 \, the second terms on 
the right-hand sides of equations (A60) and (A61) become 
negligibly small very slowly as n increases, even though 
they do become negligibly small as n approaches infinity. 
If \C 2 \ and \C 3 \ are equal, these terms never become negligi- 
bly small. Thus, the problem of circumventing this slow 
convergence or apparent lack of convergence arises. 

A satisfactory method for coping with these conditions is 
to combine linearly the results of the last two cycles of the 
series of iteration cycles that have been performed. For best 
results in an actual problem, not less than the third and 
fourth cycles should be used for this purpose in order to 
reduce as much as practicable the effects of all higher-order 
components. 

The following formulas for combining the results of the 
last two cycles are based on the assumption that the assumed 


/ rw a2 m +w g2 ( ” +I) \ __ ( rw a2 {n) Arw a 2 (n+1) \ .. . 

Vr , w a2 ( ”^ 1) +w o2 < ” > A KirWa^-v J r w a2 in) ) 2 


exists, in which stations 1 and 2 must be different or may 
be the same, depending on whether the w functions on the 
left-hand side are the same or different types of components 
than those on the right-hand side. The two values of r 
that satisfy equation (A63) are 


_^(n-l),(re+U 
r= — ,- rdz- 


in which 


2A (n_1),(re) 


J^- n ~ 1), (re) . 


(»>,(»+ 1). 


(n + 1) __ 


A (ti — 1) , (n-fl) 
7? = _ L 

2A (n_1),(71) 31 


y^4u-u,o+i>\ 

2 j^WAn+l) 

\2A^ n ~ 1,An) ) 


(w« 2 (n_1) )i 

(W o2 ( ' 1 - 1, ) 2 


(w a2 (n, )i 

(w a2 M ) 2 


(w a2 (n) )i 

(W a2 (n) ) 2 


( , ia a2 (n+1) ) 1 

(W a2 (re + 1) ) 2 


{w a2 i - n - l) ) l 

(w a2 ^) 2 


(w a2 (re+1, )i 

(w a2 ^) 2 


of R are 



1 

^2 


2 A ( - 


i ).(«> 




l),(n) 


(A 6 4) 


(A65) 
(A 6 6) 
(A67) 


(A68) 
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These values of R are equal to C 2 and C z , and the corre- 
sponding values of r, when placed in the expression 

rw a2 in} +w a2 < - n+1) (A69) 

give modes of the shapes w 2 —w x and w 3 — w u When \C 2 \ is 
nearly equal to \C Z \, the appropriate set of R and r to give the 
lower transformed mode w 2 —Wy is evident. When \C 2 \ and 
|<7 3 | are equal, the mode obtained by equation (A69) with 
either value of r may be used as the transformed second 
mode, but the trends of the eigenvalues that have been or will 
be determined at other values of the reduced frequency k 
may be used as a guide in making the selection that fits the 
trend. 

In actual computations, one further cycle of iteration be- 
ginning with an assumed mode given by expression (A69) 
should be carried out to assess the extent to which the func- 
tions w a2 in ~ 1 ' > , Wa 2 (n) , and w a2 (n+1) are free of all except the 
two components of the types appearing in equations (A60) 
and (A61). If the ratios of this cycle are not reasonably 
constant, the unwanted components still present have to be 
removed by carrying out another cycle of iteration and 
again applying equations (A64) and (A 68 ). 

The method just described is clearly applicable in the 
general cases \G n \ = \C n+l \ or \C n \ ~ \C n+l \. 

Eigenvalues having equal moduli include the special case 


of identical eigenvalues. As a basis for discussion let it be 

assumed that 


|Cil>|C 2 | = |Ci,|>|<7 4 |> . . . 

(A70) 

and that 


C 2 —C z — C 2Z 

(A71) 


The significance of the occurrence of these two identical 
eigenvalues is that the wing system may oscillate with the 
same frequency and artificial damping in any of an infinite 
number of modes, any two of which are linearly independent 
of each other and of the first, fourth, and higher modes. 
This infinite number of possible modes (all corresponding to 
C 23 ) are the infinitely many linear combinations of two 
basic linearly independent modes that are necessary and 


sufficient in combination with the first, fourth, and higher 
modes to describe an arbitrary displacement of the wing 
system. Clearly, only two linearly independent modes cor- 
responding to the double eigenvalue C 2Z are required for 
analytical purposes. These two are designated w 2 and w z as 
before but with the reservation that w 2 and w z must be 
derivable as two differing linear combinations of a single 
basic pair of linearly independent modes that also correspond 
to C 23 - 

Equations (A20) and (A 21 ) are replaced in the present case 
by 

lim w a2 (,l+1) = lim C 2Z n [b 2 (w 2 —Wi) + b z (w z —w 1 )] (A72) 

71— >co 71 — >00 

and 

y//J (« + l) 

lim (A73) 

71 — >ao lU a 2 

Equation (A27) is replaced by 
lim Wt,i (n) 

lim Wa 2 (n+1) +^j = lim C 23 n (b 2 w 2 -i- b 3 w 3 ) (A74) 

n — >00 t-'i 1 71 — ► co 

cw 

The transformed second mode (equation (A72)) is in this 
case a linear combination of the first three eigenfunctions, 
and the so-called true second mode is actually a linear 
combination of the second and third eigenfunctions. 

If the iterative transformation procedure is now applied 
in the regular way to determine the transformed third mode, 
the third eigenvalue, and the true third mode, the results 
will be as follows: The transformed third mode will be, like 
the transformed second mode, a linear combination of the 
first three eigenfunctions but will be linearly independent 
of the transformed second mode. The so-called true third 
mode will be, like the so-called true second mode, a linear 
combination of the second and third eigenfunctions and will 
be linearly independent of the so-called true second mode. 
The results will also include a second determination of the 
double eigenvalue C 23 . It may therefore be concluded that 
the iterative transformation procedure is valid and sufficient 
in all cases of eigenvalue multiplicity. 


APPENDIX B 

THE COMPLEX STIFFNESS FOR BEAMS WITH STRUCTURAL DAMPING 


The familiar concept of a complex force K(l-\~ig)s in 
simple (one-degree-of-freedom) vibrating systems having 
structural damping may be easily extended to continuous 
vibrating systems such as beams and airplane wings. The 
quantity K is the elastic-spring constant, s is the displace- 
ment, Ks is the elastic-spring force, and Kgs is the structural- 
damping force. 

For a beam in flexure, the stiffness of the fibers is given by 
the modulus of elasticity E, which is analogous to the quantity 
K for the spring. The elastic stress at any point of the cross 
section is given by eE where e is the strain which is analogous 
to the displacement s. Then the complex stress at any 


point of the cross section of a beam with structural damping 
is E(l+ig)e. The complex bending moment corresponding 
to this stress, obtained in the usual way by integration of 

cL^i] 

the moment of the stresses over the section, is EI(l -}-ig ) 

This result leads to the concept of a complex stiffness 
El (l +ig v ) for beams in flexural vibration with structural 
damping. Similarly, the complex stiffness of beams in tor- 
sional vibration with structural damping is GJ(\ -\-ig^). 
The subscripts y and <f> indicate that the structural-damping 
coefficient g may have a different value for torsional vibra- 
tions than it has for flexural vibrations. Both g v and g+ may 
be functions of the spanwise position x. 


APPENDIX C 


FORMULAS FOR EQUIVALENT CONCENTRATIONS AND INCREMENTS OF TORQUE 


The formulas used in the numerical examples for com- 
puting equivalent concentrated loads and curvatures are 
those that have been derived in references 7 and 8. For the 
concentration at an end station the formula is 

Pi=^(7pi+6p 2 —p3) (Cl) 

At an intermediate station 

(Pi-hl0p 2 +p 3 ) (C2) 

The significance of the quantities used in formulas (Cl) and 
(C2) is shown in sketch 1 . 


Distributed -load curve 



Sketch 1. 


These formulas are based on the assumption that the 
distributed-load (or curvature) curve is a series of second- 
degree parabolic arcs. When applied to distributed flexural 
loads, the formulas give concentrations which produce the 
same bending moments in the wing at all the selected stations 
as the distributed load. The formulas may be correctly applied 
to distributed torsional loads only if GJ is constant over each 
bay. In this case the formulas give concentrations which 


produce the same torsional displacement at all the selected 
stations as the distributed load. For a station placed at a 
discontinuity in ordinate or slope, formula (Cl) must be 
applied to both the left and the right of the station and the 
results added. 

The formulas for obtaining increments of area beneath a 
curve of distributed torques are derived in reference 8. 
These formulas are based as before on approximating second- 
degree parabolas. They are given here in a slightly different 
form which is better adapted to present uses. Thus 

(2i+4g 2 +23)+^ (q x — &) (C3) 

(Si + 4g 2 +(Z 3 )— £ (ffi — 2a) (C4) 


where the significance of A : and A 2 and of qi, q 2 , and q 3 is 
shown by sketch 2. 



The ordinate at a discontinuity should not be used as the 
middle one of the three ordinates selected for use in formulas 
(C3) and (C4). The formulas are valid only where the three 
ordinates are connected by a continuous curve. 
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TABLE 1. — ITERATION TO OBTAIN FIRST COUPLED MODE 
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[Common factors for each column are 
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Torsion (with GJ constant over each bay): (a) First cycle 
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FOR *:= <o FOR WING WITH CONCENTRATED MASS 
given under the column headings ] 



5 4 3 2 1 

Station 

(a) First cycle — Continued (b) Second cycle (c) Third cycle 
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TABLE 1.— ITERATION TO OBTAIN FIRST COUPLED MODE FOR k= <x> 

Torsion (with GJ variable) : (a) 


FOR WING WITH CONCENTRATED MASS— Concluded 
First cycle 
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TABLE 2.— ITERATION TO OBTAIN TRANSFORMED SECOND MODE FOR k= « FOR WING WITH CONCENTRATED MASS 

[Common factors for each column are given under the column headings] f 


Flexure: (a) First cycle 


Station 

(b) Second cycle 


(c) Third cycle 
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TABLE 4.— AERODYNAMIC-INERTIA FORCE COEFFICIENTS FOR VARIOUS VALUES OF k FOR EXAMPLE WING 

[Common factors for each column are given under the column headings] 
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A = 0.1443 FOR WING WITH CONCENTRATED MASS 
given under the column headings] 
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(a) First cycle — Continued 
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TABLE 5— ITERATION TO OBTAIN FIRST MODE FOR 


rst cycle — Continued 
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k = 0.1443 FOR WING WITH CONCENTRATED MASS— Continued 


(a) First cycle — Concluded 
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TABLE 5— ITERATION TO OBTAIN FIRST MODE FOR 


Flexure: (b) Second cycle 
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Torsion: (b) Second cycle 


(c) Third cycle 
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Xo 4 7 1 



Ely. C 

-0. 00660 

0. 00399 

- 1. 570 

1. 594 

00680 

. 00300 

-1. 690 

1. 350 

-. 00716 

. 00068 

-1. 963 

. 764 

-. 00358 

. 00036 

-. 982 

. 387 

0 

0 

0 

0 


1 

2 

3 

4 

5 

<£ir (2> 


<£ir (3) 

<t>n <3> 

+?-4>17 <3) 
<Ai72 (2> -\- i<t>U ( ' 2) 




Xo 4 7 

1 




Ely C 

-0. 00382 

0. 00466 

- 1. 407 

1. 601 

353+ 12.lt 

00407 

. 00358 

- 1. 537 

1. 353 

_ _ 

-. 00457 

. 00126 

- 1. 829 

. 773 

414—55 .It 

-. 00228 

. 00064 

-. 912 

. 392 


0 

0 

0 

0 



2 

2 




£(<W 5) +*>i/ (5) ) 


= (269.5-82.2 1) 


Xq 4 T _1 

Ely C’ 


/423000 
V 269.5 


39.7 radians per second; 
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k = 0.1443 FOR WING WITH CONCENTRATED MASS— Concluded 


(c*) Third cycle — Concluded 


(d) Fourth cycle 


5 

!/ut (i) \riyu u) 


!hn {1) + /i/i/ <3> 


Xo 4 r i 


Ely C 



267.5-82.89/ 


268-83.2/ 


(e) 'I'liird cycle — Concluded 


5 


4>in w +/4*i7 c<) 


4*177 <3) +/<#>17 (,) 


X„<7 1 


Ely C 


282-71.5 i 


282-79.8 1 


1 

2 

3 

4 

5 

6 

t/it7 <4) 

!/u u) 

2/!77 (5) 

7/l7 <5) 

20ft <5) -rR/i7 (5) 

S(7/lft (S) -f-/yi7 (5) ) 



W 4) + U/u H) 

Z(yi77 (<) + iyu M ) 

b 


Xo^iry 1 

Xo'-y 1 



Ely C 

B//. C 

1. 000 

o ; 

268. 6 


- 82. 56 

268.6-82.56/ 


. 569 

. 0030 

153. 0 


-46. 29 


j 

. 194 

. 0027 

52. 4 


- 15. 28 

270-82.5/ 

269-82.8/ 

. 054 

. 0011 

14. 7 


-4. 18 


! 

0 

0 

0 


0 


| 


(d) Fourth cycle 


1 

2 

3 

4 

5 

6 

01/7 (4) 

4*17 (4) 

4>./7 (5) 

*17<® 

4>1B (5) 'T 74*17 (5) 
4 , i/? (4) + /4>i7 (4) 

£(4>ia (5) ~r/0i7 (5) ) 
2 ( 4*177 <4) T /4>I7 <4) ) 


Xo 4 Y 1 j 

Ely C 

-0. 00704 
00720 
-. 00751 
-. 00376 
0 

0. 00379 
. 00281 
. 00052 
. 00028 

0 

— 1. 591 
-1. 711 
- 1. 977 
-. 988 
0 

1. 589 
1. 340 
. 758 
. 384 

0 

270-80.8/ 

269-82.4/ 

1 

« 

270—81.5/ 


-82.2 1 
»»«=* 269.5 =-°- 305 ; ®i=3 


39.7 

0.1443 


= 91.6 feet per second 
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TABLE 6 —ITERATION TO OBTAIN TRANSFORMED SECOND 



[Common factors for each column are 


Station 


Flexure: (a) First cycle 
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AN ITERATIVE TRANSFORMATION PROCEDURE FOR NUMERICAL SOLUTION OF FLUTTER PROBLEMS 
MODE FOR * = 0.1443 FOR WING WITH CONCENTRATED MASS 
given under the column headings] 

/[ ! © ! I 

5 4 3 2 1 


Station 


(b) Second cycle 


1 

2 

3 

4 

5 

6 

7 

8 

9 

y a iu l3) 


VbR U) 

ybi m 

VblR <2) 

Vbii (i) 

Va2R (3) 


y a 2R w + ; 

y a zR li) -r iyail™ 


b 



X 0 *by 1 



\o*y 1 





Ely C 



Ely C 

0 

0 

- 888. 5 

104. 65 

888. 5 

- 104. 65 

0 

o 


-. 188 

. 064 

-510. 8 

61. 60 

504. 5 

-57. 1 

— 6. 3 



117 

. 090 

- 177. 0 

22. 28 

172. 5 

-17. 7 

-4. 5 

4 6 

43—6 U 

— .041 

. 036 

-50. 0 

6. 45 

48. 4 

-4. 64 

-1. 6 

1. 81 


0 

0 

0 

0 

0 

0 

0 

0 



(b) Second cycle 


I 

2 

3 

4 

5 

6 

7 

S 

9 

y><i2 tt (2) 


<t>bR (2) 

<£w (2) 

$mk C2) 

4>m/ {2) 

■ 

QalR^ 

0o2/ (3) 

i 

0a?/? (3J + 6AoS/ <3} 
-T-J0a2/ <a 






Ao'y 1 
Ely C 



1 

| 

j 

1. 000 

0 

30. 80 

- 6. 92 

- 5. 88 

4. 05 

24. 92 

— 2. 87 

24.92-2.87/ j 

. 886 

— . 0153 

28. 36 

-6. 33 

-6. 12 

3. 21 

22. 24 

-3 12 

. 577 

- . 0426 

21. 60 

-4. 59 

-6. 61 

I. 24 

14. 99 

-3. 35 

26.2 -3.86/ | 

. 298 

-. 0221 

11. 01 

-2. 36 

-3. 30 

. 63 

7. 71 

— 1. 73 

0 

0 

0 

0 

0 

0 

0 

0 
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TABLE 6— ITERATION TO OBTAIN TRANSFORMED SECOND MODE 
Flexure: (c) Third cycle 




Torsion: (c) Third cycle 
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FOR k= 0.1443 FOR WING WITH CONCENTRATED MASS— Concluded 


(d) Fourth cycle 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

y a tR lt> 

*/o27 (4) 

ybR U) 

2 /w (4) 

2/mr (4) 

ybu li) 

y«2« (6) 

2/ 0 27<» 

y„ 2 R CS) +iy a 2 i < -’ j) 
2 / 0 2 r (4) +iy.u U) 

Z(yo2R <5) +M/a27 <5) ) 

Z(2/o2R (4) +*yo2J C4) ) 

6 




Xo 4 6y 1 




Xo 4 Y 1 





EIn C 




EliX C 

0 


-941. 0 

165. 58 

941. 0 

- 165. 58 

0 

0 



281 

. 184 

-543. 3 

98. 06 

534. 9 

-91. 8 

-8. 4 

6. 3 

31 — 2.0* 


252 

. 165 

- 190. 0 

35. 82 

183. 1 

- 29. 50 

-6. 9 

6. 32 

31 — 5.0* 

31-3.41* 

-. 094 


-54. 0 

10. 42 

51. 3 

-7. 93 

-2. 7 

2. 49 

31 — 3.5* 


0 

0 

0 

0 

0 

0 

0 

0 




(d) Fourth cycle 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

0a2R <4> 

0a27 <4) 

0fcfl (4> 

$ft/ (4) 

061 R (4) 

06i7 (4) 

0a2R (S) 

0 a 2 !f“> 

0a2R <S) +*0a27 <5) 
0a2R (4 * + *0a27* 4) 

£ (0a2R (5) + f0a27 (5) ) 
£ (0a2R (4) + *0a27 (4> ) 


Xo 4 y 1 
EIn C 

1. 000 

0 

33. 16 

-8. 75 

-6. 03 

4. 68 

27. 13 


-4. 07 

27.13-4.07* 


. 901 

— . 0256 

30. 70 

-8. 15 

-6. 33 

3. 80 

24. 37 


-4. 35 

27.1- 4.05* 


. 624 

-. 0749 

23. 92 

-6. 30 

-6. 98 

1. 73 

16. 94 


-4. 57 

27.7- 4.02* 

27.3-4.05* 

. 322 

-. 0384 

12. 20 

-3. 24 

-3. 49 

. 88 

8. 71 


-2. 36 

27.5- 4.06* 


0 

0 

0 

0 

0 

0 

0 


0 




cj 2 = 121.0 radians per second; g a 2 = —0.1308; *j 2 = 280 feet per second 










Cn^GOtOH-i | S’ || Cn CO tO I— ‘ 
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TABLE 7.— ITERATION TO OBTAIN TRANSFORMED THIRD 



Flexure: (a) First cycle 



Torsion: (a) First cycle 
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MODE FOR A: = 0.1443 FOR WING WITH CONCENTRATED MASS 
given under the column headings] 


I 

5 


® L 

3 2 

Station 

(b) Second cycle 


I 


1 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

2/o3 K <2) 

Vrt / t2) 

VbR (i) 

Vbi {2) 

VblR {2) 

Vbii™ 

?/6a2K <2) 

yba2I (2) 

Va 3R (W 

2/a3J t3) 

2/a3fl (3> +i2/«3i (3) 

2/o3s (8 +iy a 3r (2) 







1 




Xo 4 7 1 

L 





Bln 

C 




Eln C 

0 

0 

484. 0 

~ 83. 49 

-484. 0 

83. 49 

0 

0 

0 

0 


. 967 

013 

296. 2 

-48. 57 

-274. 8 

46. 14 

-4. 5 

2. 08 

16. 9 

-. 35 

17.5 — 0.13 i 

1. 000 

0 

113. 0 

-16. 82 

-93. 9 

14. 80 

-3. 8 

2. 23 

15. 3 

. 21 

15.3 + 0.211 

. 379 

. 004 

34. 4 

-4. 73 

-26. 3 

3. 97 

-1. 4 

. 88 

6. 6 

. 12 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 



(l>) Second cycle 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

0a3 S <2) 

0o3/ (2) 

06R <2) 

0i/ (2) 

061R l2) 

06 U (8 

<f>ba2R^ 

06a27 (2) 

0a3R <3) 

0o3/ {3) 

0O3R (3) + f0a3/ <3) 
0a 3R (2) + 10a3/ <2) 







X 0 4 7 

L 










Eln C 




0 

0 

mm 



-2. 392 

12. 87 

0. 752 

0 

0 


-. 060 

-. 0019 




- 1. 938 

11. 61 

. 358 

-1. 15 

-. 059 


-. 241 

-. 0289 


. 858 

3. 59 

-. 872 

8. 23 

-. 460 

-4. 18 

— . 474 

17. 4-0.116/ 

-. 113 

-. 0125 


. 472 

1. 80 

445 

4. 23 

-. 242 

-2. 04 

-. 215 

18. 0-0.093/ 

0 

0 

mm 

0 

0 

0 

0 

0 

0 

0 
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TABLE 7.— ITERATION TO OBTAIN TRANSFORMED THIRD MODE 


Flexure: (c) Third cycle 
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FOR *=0.1443 FOR WING WITH CONCENTRATED MASS— Concluded 

(d) Fourth cycle 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

2/o3R (4> 

y a ti w 

Vb r <4) 

Vbi w 

S/hb (<) 


yba 2 K li) 

Vb« 2I W 

y a 3s (W 

y a u (S) 

yoas m +iy a 3i w 


E(2/«3K W) + W«) 

y a m w + iy a u w 


'L{y a iR w +iy a u w ) 






Xo 4 £>7 1 





Ao 4 7 1 


b 




Ely C 





Ely C 

0 

0 

494. 8 

-86. 24 

-494. 8 

86. 24 

0 

0 

0 

0 

14. 3 + 0.46* 



. 952 

-. 033 

300. 3 

-50. 13 

-281. 6 

47. 75 

-5. 1 

2. 35 

13. 6 

-. 03 


14. 53 + 0.547* 

1. 000 

0 

115. 2 

- 17. 33 

-96. 2 

15. 31 

-4. 3 

2. 53 

14. 7 

. 51 

14. 7+0.51* 


. 390 

-. 005 

34. 4 

-4. 87 

-26. 9 

4. 12 

-1. 7 

1. 00 

5. 8 

. 25 

14. 8 + 0.83* 



0 

0 

0 

0 

0 

0 

0 

0 

0 

0 





(d) Fourth cycle 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

4>a3R (4) 

4>a31 W 

<t>bR W 

<Ai>/ <4) 

</>mr (4) 


<t>ba 2R <4) 

<t>ba21 W 

4>a3R (5) 

<£a3 f W) 

<£a3R (5) +*4>a3Z <5 ‘' 
<#>a3R (4) + ?0a37 <4) 

£ (0a3R (6> + i<t>o31 i6) ) 
X (<£a3R (4) +*4>a3/ (4) ) 


AnW 1 
ElyC 

0 

-. 083 
-. 304 
-. 150 
0 

0 

-. 0036 
-. 0326 
-. 0144 
0 

- 17. 68 
-17. 70 
- 17. 53 
-8. 86 
0 

1. 629 
1. 511 
. 843 
. 470 

0 

3. 18 
3. 33 
3. 67 
1. 84 
0 

-2. 460 
- 1. 992 
-. 901 
-. 459 
0 

14. 50 
13. 09 
9. 28 
4. 77 
0 

0. 831 
. 385 
-. 530 
-. 278 
0 

0 

-1. 28 
-4. 58 
-2. 25 
0 

0 

-. 096 
-. 588 
-. 267 
0 

15.4 + 0.49* 
15.1 + 0.32* 
15.1 + 0.33* 

15.12+0.343* 


0)3=168.9 radians per second; p a8 =0.030; w 3 = 390 feet per second 




TABLE 8.— COMPUTATION OF TRUE SECOND MODE FOR ifc=0.1443 FOR WING WITH CONCENTRATED MASS 

^Flexural functions are in terms of b ; torsional functions are in radians; l = F 12 = 8.65— 1.600i‘J 

/ 

5 




3 

Station 


Flexure 


Station 


1 (A) 


2 


3 


4 


5 



Vbia w + iybu w = 
Fuiyut+iyu) 


941.0- 165.58* 


ViR + iyn 


108.8-U0.96z | 
61.7+0.85i 
21.06+0.49i 
5.90+0.18i 
0 


W’ + iW® 


0 

■8.4+6.31 
■6.9+6.32i 
■ 2.7+2.49i 
0 


2/2 R + 12/27 


108.8 + 0.96i 
53.3 + 7.15i 
14.2 + 6.8H' 
3.2+2.67i 
0 


Torsion 
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TABLE 9.— COMPUTATION OP TRUE THIRD MODE FOR £=0.1443 FOR WING WITH CONCENTRATED MASS 

c c 

[Flexural functions are in terms of b ; torsional functions are in radians; 1 = Fi 3 = 16.97 - 6.09/; f?- 1= F 23 = 0.940- 0.313*'] 


5 4 3 2 1 

Station 


/ 

/ 

/, 


Flexure 


Station 



Station 



1 

2 

3 

l 

4 5 

6 

7 

8 

9 

10 

2/Mft (4) + «2/Mi (4) = 

F n{yujr\~iyni)-\- 

Pn{ynR-\-iyiti) 

2/6a2S t4) = 
iyba2I U) + 
Ptl(ya 2 R -\-iy a 2 l) 

P 23 (?/l2je+ i 2/12/) 

F]3(*/llK+ lynl) 2 /l 2 ft+ 0 /l 2 / 

VnR + iyui 

ym+iyu 

y a2R + /l/a27 

2/a3S (5) +iy a 3i i5) 

ViR^iy-n 

— 494.8+86.24/ 


56.3+12.30* 

-551.1 + 73.94/ 50.1 + 29.7/ 

| 

-30.2-6.5/ 

j 19.9 + 23.2/ j 
11.26+13.23/ 
3.80+4.56/ 
1.06+1.28/ 

0 

0 

-5.56 + 0.64/ 
-4.87+1.05/ 
- 1.91 + 0.42/ 
0 

0 

13.6-0.03/ 

14.7+0.51/ 

5.8+0.25/ 

0 

19.9+23.2/ 

19.3+13.84/ 

13.6+6.12/ 

5.0+1.95/ 

0 

















! 

. 



Torsion 


1 

2 

3 

4 

5 

6 7 

8 

9 

10 

<t>b\R^ + i<t>b\I^ = 

P I3(</>llB + /4>ll/) + 

Fn{4>\2R~\"><t>m) 

4>6a2fi ( ' 1) + 
/</>6 o2/ <4) = 

FiS {<t>a2R + i<t>a2l) 

F 23(<£l2R+ *012/) 

F i3(0hr + ?0h/) 

012R+?012/ 

0nR+/0n/ 0 ir + /0u 

0a2R + ?0a27 

0 o 3R <5) + /0o3J ( ® 

03R + /037 


14.50+0.831/ 





— 0.227—0.089/ 
-0.208-0.113/ 
-0.162-0.164/ 
-0.081-0.081/ 
0 

\ 13.58 + 5.43/| 
12.37+4.55/ 
9.03+2.46/ 
4.64+1.25/ 
0 

0 

-1.28-0.096/ 

-4.58-0.588/ 

-2.25-0.267/ 

0 

13.35+5.34/ 
10.88+4.34/ 
4.29+1.71/ 
2.31 + 0.90/ 
0 



| 




















i 
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